/[escript]/trunk/escript/py_src/generatediff
ViewVC logotype

Annotation of /trunk/escript/py_src/generatediff

Parent Directory Parent Directory | Revision Log Revision Log


Revision 588 - (hide annotations)
Fri Mar 10 04:45:04 2006 UTC (13 years, 7 months ago) by gross
Original Path: trunk/escript/py_src/generateutil
File size: 150465 byte(s)
1D and 3D tests for eigenvalues_and_eigenvector added.
1 jgs 154 #!/usr/bin/python
2 gross 283 # $Id$
3 jgs 154
4     """
5     program generates parts of the util.py and the test_util.py script
6     """
7 gross 157 test_header=""
8     test_header+="import unittest\n"
9     test_header+="import numarray\n"
10     test_header+="from esys.escript import *\n"
11     test_header+="from esys.finley import Rectangle\n"
12     test_header+="class Test_util2(unittest.TestCase):\n"
13     test_header+=" RES_TOL=1.e-7\n"
14     test_header+=" def setUp(self):\n"
15     test_header+=" self.__dom =Rectangle(11,11,2)\n"
16     test_header+=" self.functionspace = FunctionOnBoundary(self.__dom)\n"
17     test_tail=""
18     test_tail+="suite = unittest.TestSuite()\n"
19     test_tail+="suite.addTest(unittest.makeSuite(Test_util2))\n"
20     test_tail+="unittest.TextTestRunner(verbosity=2).run(suite)\n"
21    
22     case_set=["float","array","constData","taggedData","expandedData","Symbol"]
23     shape_set=[ (),(2,), (4,5), (6,2,2),(3,2,3,4)]
24    
25 jgs 154 t_prog=""
26 gross 157 t_prog_with_tags=""
27     t_prog_failing=""
28     u_prog=""
29 jgs 154
30 gross 157 def wherepos(arg):
31     if arg>0.:
32     return 1.
33     else:
34     return 0.
35    
36 gross 313
37 gross 157 class OPERATOR:
38     def __init__(self,nickname,rng=[-1000.,1000],test_expr="",math_expr=None,
39     numarray_expr="",symbol_expr=None,diff=None,name=""):
40     self.nickname=nickname
41     self.rng=rng
42     self.test_expr=test_expr
43     if math_expr==None:
44     self.math_expr=test_expr.replace("%a1%","arg")
45     else:
46     self.math_expr=math_expr
47     self.numarray_expr=numarray_expr
48     self.symbol_expr=symbol_expr
49     self.diff=diff
50     self.name=name
51    
52 jgs 154 import random
53     import numarray
54     import math
55     finc=1.e-6
56    
57 gross 157 def getResultCaseForBin(case0,case1):
58     if case0=="float":
59     if case1=="float":
60     case="float"
61     elif case1=="array":
62     case="array"
63     elif case1=="constData":
64     case="constData"
65     elif case1=="taggedData":
66     case="taggedData"
67     elif case1=="expandedData":
68     case="expandedData"
69     elif case1=="Symbol":
70     case="Symbol"
71     else:
72     raise ValueError,"unknown case1=%s"%case1
73     elif case0=="array":
74     if case1=="float":
75     case="array"
76     elif case1=="array":
77     case="array"
78     elif case1=="constData":
79     case="constData"
80     elif case1=="taggedData":
81     case="taggedData"
82     elif case1=="expandedData":
83     case="expandedData"
84     elif case1=="Symbol":
85     case="Symbol"
86     else:
87     raise ValueError,"unknown case1=%s"%case1
88     elif case0=="constData":
89     if case1=="float":
90     case="constData"
91     elif case1=="array":
92     case="constData"
93     elif case1=="constData":
94     case="constData"
95     elif case1=="taggedData":
96     case="taggedData"
97     elif case1=="expandedData":
98     case="expandedData"
99     elif case1=="Symbol":
100     case="Symbol"
101     else:
102     raise ValueError,"unknown case1=%s"%case1
103     elif case0=="taggedData":
104     if case1=="float":
105     case="taggedData"
106     elif case1=="array":
107     case="taggedData"
108     elif case1=="constData":
109     case="taggedData"
110     elif case1=="taggedData":
111     case="taggedData"
112     elif case1=="expandedData":
113     case="expandedData"
114     elif case1=="Symbol":
115     case="Symbol"
116     else:
117     raise ValueError,"unknown case1=%s"%case1
118     elif case0=="expandedData":
119     if case1=="float":
120     case="expandedData"
121     elif case1=="array":
122     case="expandedData"
123     elif case1=="constData":
124     case="expandedData"
125     elif case1=="taggedData":
126     case="expandedData"
127     elif case1=="expandedData":
128     case="expandedData"
129     elif case1=="Symbol":
130     case="Symbol"
131     else:
132     raise ValueError,"unknown case1=%s"%case1
133     elif case0=="Symbol":
134     if case1=="float":
135     case="Symbol"
136     elif case1=="array":
137     case="Symbol"
138     elif case1=="constData":
139     case="Symbol"
140     elif case1=="taggedData":
141     case="Symbol"
142     elif case1=="expandedData":
143     case="Symbol"
144     elif case1=="Symbol":
145     case="Symbol"
146     else:
147     raise ValueError,"unknown case1=%s"%case1
148     else:
149     raise ValueError,"unknown case0=%s"%case0
150     return case
151 jgs 154
152 gross 157
153 jgs 154 def makeArray(shape,rng):
154     l=rng[1]-rng[0]
155     out=numarray.zeros(shape,numarray.Float64)
156     if len(shape)==0:
157 gross 157 out=l*random.random()+rng[0]
158 jgs 154 elif len(shape)==1:
159     for i0 in range(shape[0]):
160     out[i0]=l*random.random()+rng[0]
161     elif len(shape)==2:
162     for i0 in range(shape[0]):
163     for i1 in range(shape[1]):
164     out[i0,i1]=l*random.random()+rng[0]
165     elif len(shape)==3:
166     for i0 in range(shape[0]):
167     for i1 in range(shape[1]):
168     for i2 in range(shape[2]):
169     out[i0,i1,i2]=l*random.random()+rng[0]
170     elif len(shape)==4:
171     for i0 in range(shape[0]):
172     for i1 in range(shape[1]):
173     for i2 in range(shape[2]):
174     for i3 in range(shape[3]):
175     out[i0,i1,i2,i3]=l*random.random()+rng[0]
176 gross 437 elif len(shape)==5:
177     for i0 in range(shape[0]):
178     for i1 in range(shape[1]):
179     for i2 in range(shape[2]):
180     for i3 in range(shape[3]):
181     for i4 in range(shape[4]):
182 gross 550 out[i0,i1,i2,i3,i4]=l*ranm.random()+rng[0]
183 jgs 154 else:
184 gross 437 raise SystemError,"rank is restricted to 5"
185 jgs 154 return out
186    
187 gross 517 def makeNumberedArray(shape,s=1.):
188     out=numarray.zeros(shape,numarray.Float64)
189     if len(shape)==0:
190     out=s*1.
191     elif len(shape)==1:
192     for i0 in range(shape[0]):
193 gross 550 out[i0]=s*int(8*random.random()+1)
194 gross 517 elif len(shape)==2:
195     for i0 in range(shape[0]):
196     for i1 in range(shape[1]):
197 gross 550 out[i0,i1]=s*int(8*random.random()+1)
198 gross 517 elif len(shape)==3:
199     for i0 in range(shape[0]):
200     for i1 in range(shape[1]):
201     for i2 in range(shape[2]):
202 gross 550 out[i0,i1,i2]=s*int(8*random.random()+1)
203 gross 517 elif len(shape)==4:
204     for i0 in range(shape[0]):
205     for i1 in range(shape[1]):
206     for i2 in range(shape[2]):
207     for i3 in range(shape[3]):
208 gross 550 out[i0,i1,i2,i3]=s*int(8*random.random()+1)
209 gross 517 else:
210     raise SystemError,"rank is restricted to 4"
211     return out
212 jgs 154
213 gross 157 def makeResult(val,test_expr):
214 jgs 154 if isinstance(val,float):
215 gross 157 out=eval(test_expr.replace("%a1%","val"))
216 jgs 154 elif len(val.shape)==0:
217 gross 157 out=eval(test_expr.replace("%a1%","val"))
218 jgs 154 elif len(val.shape)==1:
219     out=numarray.zeros(val.shape,numarray.Float64)
220     for i0 in range(val.shape[0]):
221 gross 157 out[i0]=eval(test_expr.replace("%a1%","val[i0]"))
222 jgs 154 elif len(val.shape)==2:
223     out=numarray.zeros(val.shape,numarray.Float64)
224     for i0 in range(val.shape[0]):
225     for i1 in range(val.shape[1]):
226 gross 157 out[i0,i1]=eval(test_expr.replace("%a1%","val[i0,i1]"))
227 jgs 154 elif len(val.shape)==3:
228     out=numarray.zeros(val.shape,numarray.Float64)
229     for i0 in range(val.shape[0]):
230     for i1 in range(val.shape[1]):
231     for i2 in range(val.shape[2]):
232 gross 157 out[i0,i1,i2]=eval(test_expr.replace("%a1%","val[i0,i1,i2]"))
233 jgs 154 elif len(val.shape)==4:
234     out=numarray.zeros(val.shape,numarray.Float64)
235     for i0 in range(val.shape[0]):
236     for i1 in range(val.shape[1]):
237     for i2 in range(val.shape[2]):
238     for i3 in range(val.shape[3]):
239 gross 157 out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val[i0,i1,i2,i3]"))
240 jgs 154 else:
241     raise SystemError,"rank is restricted to 4"
242     return out
243    
244 gross 157 def makeResult2(val0,val1,test_expr):
245     if isinstance(val0,float):
246     if isinstance(val1,float):
247     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
248     elif len(val1.shape)==0:
249     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
250     elif len(val1.shape)==1:
251     out=numarray.zeros(val1.shape,numarray.Float64)
252     for i0 in range(val1.shape[0]):
253     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
254     elif len(val1.shape)==2:
255     out=numarray.zeros(val1.shape,numarray.Float64)
256     for i0 in range(val1.shape[0]):
257     for i1 in range(val1.shape[1]):
258     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
259     elif len(val1.shape)==3:
260     out=numarray.zeros(val1.shape,numarray.Float64)
261     for i0 in range(val1.shape[0]):
262     for i1 in range(val1.shape[1]):
263     for i2 in range(val1.shape[2]):
264     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
265     elif len(val1.shape)==4:
266     out=numarray.zeros(val1.shape,numarray.Float64)
267     for i0 in range(val1.shape[0]):
268     for i1 in range(val1.shape[1]):
269     for i2 in range(val1.shape[2]):
270     for i3 in range(val1.shape[3]):
271     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
272     else:
273     raise SystemError,"rank of val1 is restricted to 4"
274     elif len(val0.shape)==0:
275     if isinstance(val1,float):
276     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
277     elif len(val1.shape)==0:
278     out=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1"))
279     elif len(val1.shape)==1:
280     out=numarray.zeros(val1.shape,numarray.Float64)
281     for i0 in range(val1.shape[0]):
282     out[i0]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0]"))
283     elif len(val1.shape)==2:
284     out=numarray.zeros(val1.shape,numarray.Float64)
285     for i0 in range(val1.shape[0]):
286     for i1 in range(val1.shape[1]):
287     out[i0,i1]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1]"))
288     elif len(val1.shape)==3:
289     out=numarray.zeros(val1.shape,numarray.Float64)
290     for i0 in range(val1.shape[0]):
291     for i1 in range(val1.shape[1]):
292     for i2 in range(val1.shape[2]):
293     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2]"))
294     elif len(val1.shape)==4:
295     out=numarray.zeros(val1.shape,numarray.Float64)
296     for i0 in range(val1.shape[0]):
297     for i1 in range(val1.shape[1]):
298     for i2 in range(val1.shape[2]):
299     for i3 in range(val1.shape[3]):
300     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0").replace("%a2%","val1[i0,i1,i2,i3]"))
301     else:
302     raise SystemError,"rank of val1 is restricted to 4"
303     elif len(val0.shape)==1:
304     if isinstance(val1,float):
305     out=numarray.zeros(val0.shape,numarray.Float64)
306     for i0 in range(val0.shape[0]):
307     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
308     elif len(val1.shape)==0:
309     out=numarray.zeros(val0.shape,numarray.Float64)
310     for i0 in range(val0.shape[0]):
311     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1"))
312     elif len(val1.shape)==1:
313     out=numarray.zeros(val0.shape,numarray.Float64)
314     for i0 in range(val0.shape[0]):
315     out[i0]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0]"))
316     elif len(val1.shape)==2:
317     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
318     for i0 in range(val0.shape[0]):
319     for j1 in range(val1.shape[1]):
320     out[i0,j1]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1]"))
321     elif len(val1.shape)==3:
322     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
323     for i0 in range(val0.shape[0]):
324     for j1 in range(val1.shape[1]):
325     for j2 in range(val1.shape[2]):
326     out[i0,j1,j2]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2]"))
327     elif len(val1.shape)==4:
328     out=numarray.zeros(val0.shape+val1.shape[1:],numarray.Float64)
329     for i0 in range(val0.shape[0]):
330     for j1 in range(val1.shape[1]):
331     for j2 in range(val1.shape[2]):
332     for j3 in range(val1.shape[3]):
333     out[i0,j1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0]").replace("%a2%","val1[i0,j1,j2,j3]"))
334     else:
335     raise SystemError,"rank of val1 is restricted to 4"
336     elif len(val0.shape)==2:
337     if isinstance(val1,float):
338     out=numarray.zeros(val0.shape,numarray.Float64)
339     for i0 in range(val0.shape[0]):
340     for i1 in range(val0.shape[1]):
341     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
342     elif len(val1.shape)==0:
343     out=numarray.zeros(val0.shape,numarray.Float64)
344     for i0 in range(val0.shape[0]):
345     for i1 in range(val0.shape[1]):
346     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1"))
347     elif len(val1.shape)==1:
348     out=numarray.zeros(val0.shape,numarray.Float64)
349     for i0 in range(val0.shape[0]):
350     for i1 in range(val0.shape[1]):
351     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0]"))
352     elif len(val1.shape)==2:
353     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
354     for i0 in range(val0.shape[0]):
355     for i1 in range(val0.shape[1]):
356     out[i0,i1]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1]"))
357     elif len(val1.shape)==3:
358     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
359     for i0 in range(val0.shape[0]):
360     for i1 in range(val0.shape[1]):
361     for j2 in range(val1.shape[2]):
362     out[i0,i1,j2]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2]"))
363     elif len(val1.shape)==4:
364     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
365     for i0 in range(val0.shape[0]):
366     for i1 in range(val0.shape[1]):
367     for j2 in range(val1.shape[2]):
368     for j3 in range(val1.shape[3]):
369     out[i0,i1,j2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1]").replace("%a2%","val1[i0,i1,j2,j3]"))
370     else:
371     raise SystemError,"rank of val1 is restricted to 4"
372     elif len(val0.shape)==3:
373     if isinstance(val1,float):
374     out=numarray.zeros(val0.shape,numarray.Float64)
375     for i0 in range(val0.shape[0]):
376     for i1 in range(val0.shape[1]):
377     for i2 in range(val0.shape[2]):
378     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
379     elif len(val1.shape)==0:
380     out=numarray.zeros(val0.shape,numarray.Float64)
381     for i0 in range(val0.shape[0]):
382     for i1 in range(val0.shape[1]):
383     for i2 in range(val0.shape[2]):
384     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1"))
385     elif len(val1.shape)==1:
386     out=numarray.zeros(val0.shape,numarray.Float64)
387     for i0 in range(val0.shape[0]):
388     for i1 in range(val0.shape[1]):
389     for i2 in range(val0.shape[2]):
390     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0]"))
391     elif len(val1.shape)==2:
392     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
393     for i0 in range(val0.shape[0]):
394     for i1 in range(val0.shape[1]):
395     for i2 in range(val0.shape[2]):
396     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1]"))
397     elif len(val1.shape)==3:
398     out=numarray.zeros(val0.shape,numarray.Float64)
399     for i0 in range(val0.shape[0]):
400     for i1 in range(val0.shape[1]):
401     for i2 in range(val0.shape[2]):
402     out[i0,i1,i2]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2]"))
403     elif len(val1.shape)==4:
404     out=numarray.zeros(val0.shape+val1.shape[3:],numarray.Float64)
405     for i0 in range(val0.shape[0]):
406     for i1 in range(val0.shape[1]):
407     for i2 in range(val0.shape[2]):
408     for j3 in range(val1.shape[3]):
409     out[i0,i1,i2,j3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2]").replace("%a2%","val1[i0,i1,i2,j3]"))
410     else:
411     raise SystemError,"rank of val1 is rargs[1]estricted to 4"
412     elif len(val0.shape)==4:
413     if isinstance(val1,float):
414     out=numarray.zeros(val0.shape,numarray.Float64)
415     for i0 in range(val0.shape[0]):
416     for i1 in range(val0.shape[1]):
417     for i2 in range(val0.shape[2]):
418     for i3 in range(val0.shape[3]):
419     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
420     elif len(val1.shape)==0:
421     out=numarray.zeros(val0.shape,numarray.Float64)
422     for i0 in range(val0.shape[0]):
423     for i1 in range(val0.shape[1]):
424     for i2 in range(val0.shape[2]):
425     for i3 in range(val0.shape[3]):
426     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1"))
427     elif len(val1.shape)==1:
428     out=numarray.zeros(val0.shape,numarray.Float64)
429     for i0 in range(val0.shape[0]):
430     for i1 in range(val0.shape[1]):
431     for i2 in range(val0.shape[2]):
432     for i3 in range(val0.shape[3]):
433     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0]"))
434     elif len(val1.shape)==2:
435     out=numarray.zeros(val0.shape+val1.shape[2:],numarray.Float64)
436     for i0 in range(val0.shape[0]):
437     for i1 in range(val0.shape[1]):
438     for i2 in range(val0.shape[2]):
439     for i3 in range(val0.shape[3]):
440     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1]"))
441     elif len(val1.shape)==3:
442     out=numarray.zeros(val0.shape,numarray.Float64)
443     for i0 in range(val0.shape[0]):
444     for i1 in range(val0.shape[1]):
445     for i2 in range(val0.shape[2]):
446     for i3 in range(val0.shape[3]):
447     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2]"))
448     elif len(val1.shape)==4:
449     out=numarray.zeros(val0.shape,numarray.Float64)
450     for i0 in range(val0.shape[0]):
451     for i1 in range(val0.shape[1]):
452     for i2 in range(val0.shape[2]):
453     for i3 in range(val0.shape[3]):
454     out[i0,i1,i2,i3]=eval(test_expr.replace("%a1%","val0[i0,i1,i2,i3]").replace("%a2%","val1[i0,i1,i2,i3]"))
455     else:
456     raise SystemError,"rank of val1 is restricted to 4"
457     else:
458     raise SystemError,"rank is restricted to 4"
459     return out
460 jgs 154
461 gross 157
462     def mkText(case,name,a,a1=None,use_tagging_for_expanded_data=False):
463 jgs 154 t_out=""
464 gross 157 if case=="float":
465 jgs 154 if isinstance(a,float):
466     t_out+=" %s=%s\n"%(name,a)
467 gross 291 elif a.rank==0:
468 jgs 154 t_out+=" %s=%s\n"%(name,a)
469     else:
470     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
471 gross 157 elif case=="array":
472     if isinstance(a,float):
473     t_out+=" %s=numarray.array(%s)\n"%(name,a)
474 gross 291 elif a.rank==0:
475 gross 157 t_out+=" %s=numarray.array(%s)\n"%(name,a)
476     else:
477     t_out+=" %s=numarray.array(%s)\n"%(name,a.tolist())
478 jgs 154 elif case=="constData":
479     if isinstance(a,float):
480     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
481 gross 291 elif a.rank==0:
482 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
483     else:
484     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
485     elif case=="taggedData":
486     if isinstance(a,float):
487     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
488     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
489 gross 291 elif a.rank==0:
490 jgs 154 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
491     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
492     else:
493     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
494 gross 157 t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
495 jgs 154 elif case=="expandedData":
496 gross 157 if use_tagging_for_expanded_data:
497     if isinstance(a,float):
498     t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
499     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
500 gross 291 elif a.rank==0:
501 gross 157 t_out+=" %s=Data(%s,self.functionspace)\n"%(name,a)
502     t_out+=" %s.setTaggedValue(1,%s)\n"%(name,a1)
503     else:
504     t_out+=" %s=Data(numarray.array(%s),self.functionspace)\n"%(name,a.tolist())
505     t_out+=" %s.setTaggedValue(1,numarray.array(%s))\n"%(name,a1.tolist())
506     t_out+=" %s.expand()\n"%name
507     else:
508 gross 550 t_out+=" msk_%s=whereZero(self.functionspace.getX()[0],1.e-8)\n"%name
509 gross 157 if isinstance(a,float):
510 gross 550 t_out+=" %s=(1.-msk_%s)*(%s)+msk_%s*(%s)\n"%(name,name,a,name,a1)
511 gross 291 elif a.rank==0:
512 gross 157 t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a,name,a1)
513     else:
514     t_out+=" %s=msk_%s*numarray.array(%s)+(1.-msk_%s)*numarray.array(%s)\n"%(name,name,a.tolist(),name,a1.tolist())
515 jgs 154 elif case=="Symbol":
516     if isinstance(a,float):
517     t_out+=" %s=Symbol(shape=())\n"%(name)
518 gross 291 elif a.rank==0:
519 jgs 154 t_out+=" %s=Symbol(shape=())\n"%(name)
520     else:
521     t_out+=" %s=Symbol(shape=%s)\n"%(name,str(a.shape))
522    
523     return t_out
524    
525 gross 587 def mkTypeAndShapeTest(case,sh,argstr,name=""):
526 gross 157 text=""
527     if case=="float":
528 gross 587 text+=" self.failUnless(isinstance(%s,float),\"wrong type of result%s.\")\n"%(argstr,name)
529 gross 157 elif case=="array":
530 gross 587 text+=" self.failUnless(isinstance(%s,numarray.NumArray),\"wrong type of result%s.\")\n"%(argstr,name)
531     text+=" self.failUnlessEqual(%s.shape,%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
532 gross 157 elif case in ["constData","taggedData","expandedData"]:
533 gross 587 text+=" self.failUnless(isinstance(%s,Data),\"wrong type of result%s.\")\n"%(argstr,name)
534     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
535 gross 157 elif case=="Symbol":
536 gross 587 text+=" self.failUnless(isinstance(%s,Symbol),\"wrong type of result%s.\")\n"%(argstr,name)
537     text+=" self.failUnlessEqual(%s.getShape(),%s,\"wrong shape of result%s.\")\n"%(argstr,str(sh),name)
538 gross 157 return text
539 jgs 154
540 gross 157 def mkCode(txt,args=[],intend=""):
541     s=txt.split("\n")
542     if len(s)>1:
543     out=""
544     for l in s:
545     out+=intend+l+"\n"
546     else:
547     out="%sreturn %s\n"%(intend,txt)
548     c=1
549     for r in args:
550     out=out.replace("%%a%s%%"%c,r)
551     return out
552 gross 587
553    
554 gross 530 #=======================================================================================================
555 gross 587 # eigenvalues and eigen vectors 2D:
556     #=======================================================================================================
557     alltests= \
558     [ ("case0",[[0.0, 0.0], [0.0, 0.0]],[0.0, 0.0]) \
559     , ("case3",[[-1.0, 0.0], [0.0, -1.0]],[-1.0, -1.0]) \
560     , ("case5",[[-0.99999999999999967, -6.4606252205695602e-16], [-6.4606252205695602e-16, -0.99999999999999967]],[-1.0, -1.0]) \
561     , ("case6",[[0.0, 0.0], [0.0, 0.0001]],[0.0, 0.0001]) \
562     , ("case7",[[0.0001, 0.0], [0.0, 0.0]],[0.0, 0.0001]) \
563     , ("case8",[[6.0598371831785722e-06, 2.3859213977648625e-05], [2.3859213977648629e-05, 9.3940162816821425e-05]],[0.0, 0.0001]) \
564     , ("case9",[[1.0, 0.0], [0.0, 2.0]],[1.0, 2.0]) \
565     , ("case10",[[2.0, 0.0], [0.0, 1.0]],[1.0, 2.0]) \
566     , ("case11",[[1.0605983718317855, 0.23859213977648688], [0.23859213977648688, 1.9394016281682138]],[1.0, 2.0]) \
567     , ("case12",[[1.0, 0.0], [0.0, 1000000.0]],[1.0, 1000000.0]) \
568     , ("case13",[[1000000.0, 0.0], [0.0, 1.0]],[1.0, 1000000.0]) \
569     , ("case14",[[60599.311233413886, 238591.90118434647], [238591.90118434647, 939401.68876658613]],[1.0, 1000000.0]) \
570     ]
571     dim=2
572 gross 588
573    
574     alltests= \
575     [ ("case0",[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],[0.0, 0.0, 0.0]) \
576     , ("case5",[[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]],[10.0, 10.0, 10.0]) \
577     , ("case10",[[0.9, 0.0, 0.0], [0.0, 0.9, 0.0], [0.0, 0.0, 1.0]],[0.9, 0.9, 1.0]) \
578     , ("case11",[[0.9, 0.0, 0.0], [0.0, 0.97060899725040983, -0.045555123008643325], [0.0, -0.045555123008643339, 0.92939100274959041]],[0.9, 0.9, 1.0]) \
579     , ("case12",[[0.92694799760252555, 0.0, 0.044368966468320177], [0.0, 0.9, 0.0], [0.044368966468320184, 0.0, 0.97305200239747425]],[0.9, 0.9, 1.0]) \
580     , ("case13",[[1.0, 0.0, 0.0], [0.0, 0.9, 0.], [0.0, 0., 0.9]],[0.9, 0.9, 1.0]) \
581     , ("case14",[[0.92379770619813639, 0.041031106298491521, -0.011396846732439278], [0.041031106298491535, 0.97074428392640366, -0.019650012730342326], [-0.011396846732439236, -0.019650012730342337, 0.90545800987545966]],[0.9, 0.9, 1.0]) \
582     , ("case15",[[1.0, 0.0, 0.0], [0.0, 1.1, 0.0], [0.0, 0.0, 1.1]],[1.0, 1.1, 1.1]) \
583     , ("case17",[[1.0269479976025255, 0.0, 0.044368966468320309], [0.0, 1.1, 0.0], [0.044368966468320295, 0.0, 1.0730520023974743]],[1.0, 1.1, 1.1]) \
584     , ("case18",[[1.1, 0.0, 0.0], [0.0, 1.0153410887977139, -0.036038311201720394], [0.0, -0.036038311201720373, 1.084658911202286]],[1.0, 1.1, 1.1]) \
585     , ("case19",[[1.035487967756175, 0.026317079185831614, -0.039960133424212368], [0.026317079185831618, 1.0892641940924184, 0.016301362071911414], [-0.039960133424212355, 0.016301362071911431, 1.0752478381514063]],[1.0, 1.1, 1.1]) \
586     , ("case20",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]],[1.0, 2.0, 3.0]) \
587     , ("case21",[[1.0, 0.0, 0.0], [0.0, 2.7060899725040968, -0.45555123008643206], [0.0, -0.45555123008643228, 2.2939100274959037]],[1.0, 2.0, 3.0]) \
588     , ("case22",[[1.5389599520505153, 0.0, 0.88737932936638753], [0.0, 2.0, 0.0], [0.88737932936638753, 0.0, 2.4610400479494858]],[1.0, 2.0, 3.0]) \
589     , ("case23",[[3.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 3.0]) \
590     , ("case24",[[1.5928567395431172, 0.67348185484323142, -0.51356980156651744], [0.67348185484323153, 2.6000847801882254, -0.033486506584313548], [-0.51356980156651744, -0.033486506584313541, 1.8070584802686565]],[1.0, 2.0, 3.0]) \
591     , ("case25",[[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 30000.0]],[1.0, 2.0, 30000.0]) \
592     , ("case26",[[1.0, 0.0, 0.0], [0.0, 21183.286995177881, -13665.625800132779], [0.0, -13665.625800132779, 8818.7130048221279]],[1.0, 2.0, 30000.0]) \
593     , ("case27",[[8085.1298007817086, 0.0, 13310.246250831115], [0.0, 2.0, 0.0], [13310.246250831115, 0.0, 21915.870199218316]],[1.0, 2.0, 30000.0]) \
594     , ("case28",[[30000.0, 0.0, 0.0], [0.0, 1.153410887977139, -0.36038311201720391], [0.0, -0.36038311201720391, 1.8465891120228608]],[1.0, 2.0, 30000.0]) \
595     , ("case29",[[7140.1907849945546, 12308.774438213351, -3419.2256841313947], [12308.774438213351, 21223.762934183575, -5894.4478052274408], [-3419.2256841313947, -5894.4478052274408, 1639.0462808218595]],[1.0, 2.0, 30000.0]) \
596     ]
597    
598     dim=3
599    
600     alltests= \
601     [ ("case0",[[0.0]],[0.0]) \
602     , ("case1",[[1.0]],[1.0]) \
603     ]
604     dim=1
605    
606 gross 587 for case in ["constData","taggedData","expandedData"]:
607     n=0
608     while n<len(alltests):
609     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
610     tname="test_eigenvalues_and_eigenvectors_%s_dim%s_%s"%(case,dim,alltests[n][0])
611     text+=" def %s(self):\n"%tname
612     a_0=numarray.array(alltests[n][1],numarray.Float64)
613     ev_0=numarray.array(alltests[n][2],numarray.Float64)
614     if case in ["taggedData", "expandedData"]:
615 gross 588 if n+1<len(alltests):
616     a1_0=numarray.array(alltests[n+1][1],numarray.Float64)
617     ev1_0=numarray.array(alltests[n+1][2],numarray.Float64)
618     else:
619     a1_0=numarray.array(alltests[0][1],numarray.Float64)
620     ev1_0=numarray.array(alltests[0][2],numarray.Float64)
621 gross 587 n+=2
622     else:
623     a1_0=a_0
624     ev1_0=ev_0
625     n+=1
626     text+=mkText(case,"arg",a_0,a1_0)
627     text+=" res=eigenvalues_and_eigenvectors(arg)\n"
628     text+=mkText(case,"ref_ev",ev_0,ev1_0)
629     text+=mkTypeAndShapeTest(case,(dim,),"res[0]"," for eigenvalues")
630     text+=mkTypeAndShapeTest(case,(dim,dim),"res[1]"," for eigenvectors")
631     text+=" self.failUnless(Lsup(res[0]-ref_ev)<=self.RES_TOL*Lsup(ref_ev),\"wrong eigenvalues\")\n"
632     for i in range(dim):
633     text+=" self.failUnless(Lsup(matrixmult(arg,res[1][:,%s])-res[0][%s]*res[1][:,%s])<=self.RES_TOL*Lsup(res[0]),\"wrong eigenvector %s\")\n"%(i,i,i,i)
634     text+=" self.failUnless(Lsup(matrixmult(transpose(res[1]),res[1])-kronecker(%s))<=self.RES_TOL,\"eigenvectors are not orthonormal\")\n"%dim
635     if case == "taggedData" :
636     t_prog_with_tags+=text
637     else:
638     t_prog+=text
639     print test_header
640     print t_prog
641     print t_prog_with_tags
642     print test_tail
643     1/0
644     #=======================================================================================================
645 gross 550 # get slices
646 gross 536 #=======================================================================================================
647     from esys.escript import *
648 gross 550 for case0 in ["constData","taggedData","expandedData"]:
649     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,3,3)]:
650     # get perm:
651     if len(sh0)==2:
652     check=[[1,0]]
653     elif len(sh0)==3:
654     check=[[1,0,2],
655     [1,2,0],
656     [2,1,0],
657     [2,0,2],
658     [0,2,1]]
659     elif len(sh0)==4:
660     check=[[0,1,3,2],
661     [0,2,1,3],
662     [0,2,3,1],
663     [0,3,2,1],
664     [0,3,1,2] ,
665     [1,0,2,3],
666     [1,0,3,2],
667     [1,2,0,3],
668     [1,2,3,0],
669     [1,3,2,0],
670     [1,3,0,2],
671     [2,0,1,3],
672     [2,0,3,1],
673     [2,1,0,3],
674     [2,1,3,0],
675     [2,3,1,0],
676     [2,3,0,1],
677     [3,0,1,2],
678     [3,0,2,1],
679     [3,1,0,2],
680     [3,1,2,0],
681     [3,2,1,0],
682     [3,2,0,1]]
683     else:
684     check=[]
685    
686     # create the test cases:
687     processed=[]
688     l=["R","U","L","P","C","N"]
689     c=[""]
690     for i in range(len(sh0)):
691     tmp=[]
692     for ci in c:
693     tmp+=[ci+li for li in l]
694     c=tmp
695     # SHUFFLE
696     c2=[]
697     while len(c)>0:
698     i=int(random.random()*len(c))
699     c2.append(c[i])
700     del c[i]
701     c=c2
702     for ci in c:
703     t=""
704     sh=()
705     sl=()
706     for i in range(len(ci)):
707     if ci[i]=="R":
708     s="%s:%s"%(1,sh0[i]-1)
709     sl=sl+(slice(1,sh0[i]-1),)
710     sh=sh+(sh0[i]-2,)
711     if ci[i]=="U":
712     s=":%s"%(sh0[i]-1)
713     sh=sh+(sh0[i]-1,)
714     sl=sl+(slice(0,sh0[i]-1),)
715     if ci[i]=="L":
716     s="2:"
717     sh=sh+(sh0[i]-2,)
718     sl=sl+(slice(2,sh0[i]),)
719     if ci[i]=="P":
720     s="%s"%(int(sh0[i]/2))
721     sl=sl+(int(sh0[i]/2),)
722     if ci[i]=="C":
723     s=":"
724     sh=sh+(sh0[i],)
725     sl=sl+(slice(0,sh0[i]),)
726     if ci[i]=="N":
727     s=""
728     sh=sh+(sh0[i],)
729     if len(s)>0:
730     if not t=="": t+=","
731     t+=s
732     if len(sl)==1: sl=sl[0]
733     N_found=False
734     noN_found=False
735     process=len(t)>0
736     for i in ci:
737     if i=="N":
738     if not noN_found and N_found: process=False
739     N_found=True
740     else:
741     if N_found: process=False
742     noNfound=True
743     # is there a similar one processed allready
744     if process and ci.find("N")==-1:
745     for ci2 in processed:
746     for chi in check:
747     is_perm=True
748     for i in range(len(chi)):
749     if not ci[i]==ci2[chi[i]]: is_perm=False
750     if is_perm: process=False
751     # if not process: print ci," rejected"
752     if process:
753     processed.append(ci)
754     for case1 in ["array","constData","taggedData","expandedData"]:
755     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
756     tname="test_setslice_%s_rank%s_%s_%s"%(case0,len(sh0),case1,ci)
757     text+=" def %s(self):\n"%tname
758     a_0=makeNumberedArray(sh0)
759     if case0 in ["taggedData", "expandedData"]:
760     a1_0=makeNumberedArray(sh0)
761     else:
762     a1_0=a_0*1.
763    
764     a_1=makeNumberedArray(sh)
765     if case1 in ["taggedData", "expandedData"]:
766     a1_1=makeNumberedArray(sh)
767     else:
768     a1_1=a_1*1.
769    
770     text+=mkText(case0,"arg",a_0,a1_0)
771     text+=mkText(case1,"val",a_1,a1_1)
772     text+=" arg[%s]=val\n"%t
773     a_0.__setitem__(sl,a_1)
774     a1_0.__setitem__(sl,a1_1)
775     if Lsup(a_0-a1_0)==0.:
776     text+=mkText("constData","ref",a_0,a1_0)
777     else:
778     text+=mkText("expandedData","ref",a_0,a1_0)
779     text+=" self.failUnless(Lsup(arg-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"
780    
781     if case0 == "taggedData" or case1 == "taggedData":
782     t_prog_with_tags+=text
783     else:
784     t_prog+=text
785    
786     print test_header
787     # print t_prog
788     print t_prog_with_tags
789     print test_tail
790     1/0
791    
792     #=======================================================================================================
793     # (non)symmetric part
794     #=======================================================================================================
795     from esys.escript import *
796 gross 536 for name in ["symmetric", "nonsymmetric"]:
797     f=1.
798     if name=="nonsymmetric": f=-1
799     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
800     for sh0 in [ (3,3), (2,3,2,3)]:
801     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
802     tname="test_%s_%s_rank%s"%(name,case0,len(sh0))
803     text+=" def %s(self):\n"%tname
804     a_0=makeNumberedArray(sh0,s=1.)
805     r_0=(a_0+f*transpose(a_0))/2.
806     if case0 in ["taggedData", "expandedData"]:
807     a1_0=makeNumberedArray(sh0,s=-1.)
808     r1_0=(a1_0+f*transpose(a1_0))/2.
809     else:
810     a1_0=a_0
811     r1_0=r_0
812     text+=mkText(case0,"arg",a_0,a1_0)
813     text+=" res=%s(arg)\n"%name
814     if case0=="Symbol":
815     text+=mkText("array","s",a_0,a1_0)
816     text+=" sub=res.substitute({arg:s})\n"
817     res="sub"
818     text+=mkText("array","ref",r_0,r1_0)
819     else:
820     res="res"
821     text+=mkText(case0,"ref",r_0,r1_0)
822     text+=mkTypeAndShapeTest(case0,sh0,"res")
823     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
824    
825     if case0 == "taggedData" :
826     t_prog_with_tags+=text
827     else:
828     t_prog+=text
829     print test_header
830     print t_prog
831     # print t_prog_with_tags
832     print test_tail
833     1/0
834    
835     #=======================================================================================================
836 gross 530 # eigenvalues
837     #=======================================================================================================
838     import numarray.linear_algebra
839     name="eigenvalues"
840     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
841     for sh0 in [ (1,1), (2,2), (3,3)]:
842     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
843     tname="test_%s_%s_dim%s"%(name,case0,sh0[0])
844     text+=" def %s(self):\n"%tname
845     a_0=makeArray(sh0,[-1.,1])
846     a_0=(a_0+numarray.transpose(a_0))/2.
847     ev=numarray.linear_algebra.eigenvalues(a_0)
848     ev.sort()
849     if case0 in ["taggedData", "expandedData"]:
850     a1_0=makeArray(sh0,[-1.,1])
851     a1_0=(a1_0+numarray.transpose(a1_0))/2.
852     ev1=numarray.linear_algebra.eigenvalues(a1_0)
853     ev1.sort()
854     else:
855     a1_0=a_0
856     ev1=ev
857     text+=mkText(case0,"arg",a_0,a1_0)
858     text+=" res=%s(arg)\n"%name
859     if case0=="Symbol":
860     text+=mkText("array","s",a_0,a1_0)
861     text+=" sub=res.substitute({arg:s})\n"
862     res="sub"
863     text+=mkText("array","ref",ev,ev1)
864     else:
865     res="res"
866     text+=mkText(case0,"ref",ev,ev1)
867     text+=mkTypeAndShapeTest(case0,(sh0[0],),"res")
868     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
869    
870     if case0 == "taggedData" :
871     t_prog_with_tags+=text
872     else:
873     t_prog+=text
874     print test_header
875     # print t_prog
876     print t_prog_with_tags
877     print test_tail
878     1/0
879 jgs 154
880 gross 517 #=======================================================================================================
881 gross 550 # get slices
882 gross 517 #=======================================================================================================
883     for case0 in ["constData","taggedData","expandedData","Symbol"]:
884     for sh0 in [ (3,), (3,4), (3,4,3) ,(4,3,5,3)]:
885     # get perm:
886     if len(sh0)==2:
887     check=[[1,0]]
888     elif len(sh0)==3:
889     check=[[1,0,2],
890     [1,2,0],
891     [2,1,0],
892     [2,0,2],
893     [0,2,1]]
894     elif len(sh0)==4:
895     check=[[0,1,3,2],
896     [0,2,1,3],
897     [0,2,3,1],
898     [0,3,2,1],
899     [0,3,1,2] ,
900     [1,0,2,3],
901     [1,0,3,2],
902     [1,2,0,3],
903     [1,2,3,0],
904     [1,3,2,0],
905     [1,3,0,2],
906     [2,0,1,3],
907     [2,0,3,1],
908     [2,1,0,3],
909     [2,1,3,0],
910     [2,3,1,0],
911     [2,3,0,1],
912     [3,0,1,2],
913     [3,0,2,1],
914     [3,1,0,2],
915     [3,1,2,0],
916     [3,2,1,0],
917     [3,2,0,1]]
918     else:
919     check=[]
920    
921     # create the test cases:
922     processed=[]
923     l=["R","U","L","P","C","N"]
924     c=[""]
925     for i in range(len(sh0)):
926     tmp=[]
927     for ci in c:
928     tmp+=[ci+li for li in l]
929     c=tmp
930     # SHUFFLE
931     c2=[]
932     while len(c)>0:
933     i=int(random.random()*len(c))
934     c2.append(c[i])
935     del c[i]
936     c=c2
937     for ci in c:
938     t=""
939     sh=()
940     for i in range(len(ci)):
941     if ci[i]=="R":
942     s="%s:%s"%(1,sh0[i]-1)
943     sh=sh+(sh0[i]-2,)
944     if ci[i]=="U":
945     s=":%s"%(sh0[i]-1)
946     sh=sh+(sh0[i]-1,)
947     if ci[i]=="L":
948     s="2:"
949     sh=sh+(sh0[i]-2,)
950     if ci[i]=="P":
951     s="%s"%(int(sh0[i]/2))
952     if ci[i]=="C":
953     s=":"
954     sh=sh+(sh0[i],)
955     if ci[i]=="N":
956     s=""
957     sh=sh+(sh0[i],)
958     if len(s)>0:
959     if not t=="": t+=","
960     t+=s
961     N_found=False
962     noN_found=False
963     process=len(t)>0
964     for i in ci:
965     if i=="N":
966     if not noN_found and N_found: process=False
967     N_found=True
968     else:
969     if N_found: process=False
970     noNfound=True
971     # is there a similar one processed allready
972     if process and ci.find("N")==-1:
973     for ci2 in processed:
974     for chi in check:
975     is_perm=True
976     for i in range(len(chi)):
977     if not ci[i]==ci2[chi[i]]: is_perm=False
978     if is_perm: process=False
979     # if not process: print ci," rejected"
980     if process:
981     processed.append(ci)
982     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
983     tname="test_getslice_%s_rank%s_%s"%(case0,len(sh0),ci)
984     text+=" def %s(self):\n"%tname
985     a_0=makeNumberedArray(sh0,s=1)
986     if case0 in ["taggedData", "expandedData"]:
987     a1_0=makeNumberedArray(sh0,s=-1.)
988     else:
989     a1_0=a_0
990     r=eval("a_0[%s]"%t)
991     r1=eval("a1_0[%s]"%t)
992     text+=mkText(case0,"arg",a_0,a1_0)
993     text+=" res=arg[%s]\n"%t
994     if case0=="Symbol":
995     text+=mkText("array","s",a_0,a1_0)
996     text+=" sub=res.substitute({arg:s})\n"
997     res="sub"
998     text+=mkText("array","ref",r,r1)
999     else:
1000     res="res"
1001     text+=mkText(case0,"ref",r,r1)
1002     text+=mkTypeAndShapeTest(case0,sh,"res")
1003     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1004    
1005     if case0 == "taggedData" :
1006     t_prog_with_tags+=text
1007     else:
1008     t_prog+=text
1009    
1010     print test_header
1011     # print t_prog
1012     print t_prog_with_tags
1013     print test_tail
1014     1/0
1015     #============================================================================================
1016 gross 291 def innerTEST(arg0,arg1):
1017     if isinstance(arg0,float):
1018     out=numarray.array(arg0*arg1)
1019     else:
1020     out=(arg0*arg1).sum()
1021     return out
1022    
1023     def outerTEST(arg0,arg1):
1024     if isinstance(arg0,float):
1025     out=numarray.array(arg0*arg1)
1026     elif isinstance(arg1,float):
1027     out=numarray.array(arg0*arg1)
1028     else:
1029     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1030     return out
1031    
1032     def tensorProductTest(arg0,arg1,sh_s):
1033     if isinstance(arg0,float):
1034     out=numarray.array(arg0*arg1)
1035     elif isinstance(arg1,float):
1036     out=numarray.array(arg0*arg1)
1037     elif len(sh_s)==0:
1038     out=numarray.outerproduct(arg0,arg1).resize(arg0.shape+arg1.shape)
1039     else:
1040     l=len(sh_s)
1041     sh0=arg0.shape[:arg0.rank-l]
1042     sh1=arg1.shape[l:]
1043     ls,l0,l1=1,1,1
1044     for i in sh_s: ls*=i
1045     for i in sh0: l0*=i
1046     for i in sh1: l1*=i
1047     out1=numarray.outerproduct(arg0,arg1).resize((l0,ls,ls,l1))
1048     out2=numarray.zeros((l0,l1),numarray.Float)
1049     for i0 in range(l0):
1050     for i1 in range(l1):
1051     for i in range(ls): out2[i0,i1]+=out1[i0,i,i,i1]
1052     out=out2.resize(sh0+sh1)
1053     return out
1054    
1055     def testMatrixMult(arg0,arg1,sh_s):
1056     return numarray.matrixmultiply(arg0,arg1)
1057    
1058    
1059     def testTensorMult(arg0,arg1,sh_s):
1060     if len(arg0)==2:
1061     return numarray.matrixmultiply(arg0,arg1)
1062     else:
1063     if arg1.rank==4:
1064     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2],arg1.shape[3]),numarray.Float)
1065     for i0 in range(arg0.shape[0]):
1066     for i1 in range(arg0.shape[1]):
1067     for i2 in range(arg0.shape[2]):
1068     for i3 in range(arg0.shape[3]):
1069     for j2 in range(arg1.shape[2]):
1070     for j3 in range(arg1.shape[3]):
1071     out[i0,i1,j2,j3]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2,j3]
1072     elif arg1.rank==3:
1073     out=numarray.zeros((arg0.shape[0],arg0.shape[1],arg1.shape[2]),numarray.Float)
1074     for i0 in range(arg0.shape[0]):
1075     for i1 in range(arg0.shape[1]):
1076     for i2 in range(arg0.shape[2]):
1077     for i3 in range(arg0.shape[3]):
1078     for j2 in range(arg1.shape[2]):
1079     out[i0,i1,j2]+=arg0[i0,i1,i2,i3]*arg1[i2,i3,j2]
1080     elif arg1.rank==2:
1081     out=numarray.zeros((arg0.shape[0],arg0.shape[1]),numarray.Float)
1082     for i0 in range(arg0.shape[0]):
1083     for i1 in range(arg0.shape[1]):
1084     for i2 in range(arg0.shape[2]):
1085     for i3 in range(arg0.shape[3]):
1086     out[i0,i1]+=arg0[i0,i1,i2,i3]*arg1[i2,i3]
1087     return out
1088 gross 313
1089     def testReduce(arg0,init_val,test_expr,post_expr):
1090     out=init_val
1091     if isinstance(arg0,float):
1092     out=eval(test_expr.replace("%a1%","arg0"))
1093     elif arg0.rank==0:
1094     out=eval(test_expr.replace("%a1%","arg0"))
1095     elif arg0.rank==1:
1096     for i0 in range(arg0.shape[0]):
1097     out=eval(test_expr.replace("%a1%","arg0[i0]"))
1098     elif arg0.rank==2:
1099     for i0 in range(arg0.shape[0]):
1100     for i1 in range(arg0.shape[1]):
1101     out=eval(test_expr.replace("%a1%","arg0[i0,i1]"))
1102     elif arg0.rank==3:
1103     for i0 in range(arg0.shape[0]):
1104     for i1 in range(arg0.shape[1]):
1105     for i2 in range(arg0.shape[2]):
1106     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2]"))
1107     elif arg0.rank==4:
1108     for i0 in range(arg0.shape[0]):
1109     for i1 in range(arg0.shape[1]):
1110     for i2 in range(arg0.shape[2]):
1111     for i3 in range(arg0.shape[3]):
1112     out=eval(test_expr.replace("%a1%","arg0[i0,i1,i2,i3]"))
1113     return eval(post_expr)
1114 gross 396
1115     def clipTEST(arg0,mn,mx):
1116     if isinstance(arg0,float):
1117     return max(min(arg0,mx),mn)
1118     out=numarray.zeros(arg0.shape,numarray.Float64)
1119     if arg0.rank==1:
1120     for i0 in range(arg0.shape[0]):
1121     out[i0]=max(min(arg0[i0],mx),mn)
1122     elif arg0.rank==2:
1123     for i0 in range(arg0.shape[0]):
1124     for i1 in range(arg0.shape[1]):
1125     out[i0,i1]=max(min(arg0[i0,i1],mx),mn)
1126     elif arg0.rank==3:
1127     for i0 in range(arg0.shape[0]):
1128     for i1 in range(arg0.shape[1]):
1129     for i2 in range(arg0.shape[2]):
1130     out[i0,i1,i2]=max(min(arg0[i0,i1,i2],mx),mn)
1131     elif arg0.rank==4:
1132     for i0 in range(arg0.shape[0]):
1133     for i1 in range(arg0.shape[1]):
1134     for i2 in range(arg0.shape[2]):
1135     for i3 in range(arg0.shape[3]):
1136     out[i0,i1,i2,i3]=max(min(arg0[i0,i1,i2,i3],mx),mn)
1137     return out
1138     def minimumTEST(arg0,arg1):
1139     if isinstance(arg0,float):
1140     if isinstance(arg1,float):
1141     if arg0>arg1:
1142     return arg1
1143     else:
1144     return arg0
1145     else:
1146     arg0=numarray.ones(arg1.shape)*arg0
1147     else:
1148     if isinstance(arg1,float):
1149     arg1=numarray.ones(arg0.shape)*arg1
1150     out=numarray.zeros(arg0.shape,numarray.Float64)
1151     if arg0.rank==0:
1152     if arg0>arg1:
1153     out=arg1
1154     else:
1155     out=arg0
1156     elif arg0.rank==1:
1157     for i0 in range(arg0.shape[0]):
1158     if arg0[i0]>arg1[i0]:
1159     out[i0]=arg1[i0]
1160     else:
1161     out[i0]=arg0[i0]
1162     elif arg0.rank==2:
1163     for i0 in range(arg0.shape[0]):
1164     for i1 in range(arg0.shape[1]):
1165     if arg0[i0,i1]>arg1[i0,i1]:
1166     out[i0,i1]=arg1[i0,i1]
1167     else:
1168     out[i0,i1]=arg0[i0,i1]
1169     elif arg0.rank==3:
1170     for i0 in range(arg0.shape[0]):
1171     for i1 in range(arg0.shape[1]):
1172     for i2 in range(arg0.shape[2]):
1173     if arg0[i0,i1,i2]>arg1[i0,i1,i2]:
1174     out[i0,i1,i2]=arg1[i0,i1,i2]
1175     else:
1176     out[i0,i1,i2]=arg0[i0,i1,i2]
1177     elif arg0.rank==4:
1178     for i0 in range(arg0.shape[0]):
1179     for i1 in range(arg0.shape[1]):
1180     for i2 in range(arg0.shape[2]):
1181     for i3 in range(arg0.shape[3]):
1182     if arg0[i0,i1,i2,i3]>arg1[i0,i1,i2,i3]:
1183     out[i0,i1,i2,i3]=arg1[i0,i1,i2,i3]
1184     else:
1185     out[i0,i1,i2,i3]=arg0[i0,i1,i2,i3]
1186     return out
1187 gross 443
1188 gross 439 def unrollLoops(a,b,o,arg,tap="",x="x"):
1189 gross 437 out=""
1190     if a.rank==1:
1191     z=""
1192     for i99 in range(a.shape[0]):
1193     if not z=="": z+="+"
1194     if o=="1":
1195 gross 439 z+="(%s)*%s[%s]"%(a[i99]+b[i99],x,i99)
1196 gross 437 else:
1197 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i99],x,i99,b[i99],x,i99)
1198 gross 437
1199     out+=tap+"%s=%s\n"%(arg,z)
1200    
1201     elif a.rank==2:
1202     for i0 in range(a.shape[0]):
1203     z=""
1204     for i99 in range(a.shape[1]):
1205     if not z=="": z+="+"
1206     if o=="1":
1207     z+="(%s)*x[%s]"%(a[i0,i99]+b[i0,i99],i99)
1208     else:
1209 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i99],x,i99,b[i0,i99],x,i99)
1210 gross 437
1211     out+=tap+"%s[%s]=%s\n"%(arg,i0,z)
1212     elif a.rank==3:
1213     for i0 in range(a.shape[0]):
1214     for i1 in range(a.shape[1]):
1215     z=""
1216     for i99 in range(a.shape[2]):
1217     if not z=="": z+="+"
1218     if o=="1":
1219 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i99]+b[i0,i1,i99],x,i99)
1220 gross 437 else:
1221 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i99],x,i99,b[i0,i1,i99],x,i99)
1222 gross 437
1223     out+=tap+"%s[%s,%s]=%s\n"%(arg,i0,i1,z)
1224     elif a.rank==4:
1225     for i0 in range(a.shape[0]):
1226     for i1 in range(a.shape[1]):
1227     for i2 in range(a.shape[2]):
1228     z=""
1229     for i99 in range(a.shape[3]):
1230     if not z=="": z+="+"
1231     if o=="1":
1232 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i99]+b[i0,i1,i2,i99],x,i99)
1233 gross 437 else:
1234 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i99],x,i99,b[i0,i1,i2,i99],x,i99)
1235 gross 437
1236     out+=tap+"%s[%s,%s,%s]=%s\n"%(arg,i0,i1,i2,z)
1237     elif a.rank==5:
1238     for i0 in range(a.shape[0]):
1239     for i1 in range(a.shape[1]):
1240     for i2 in range(a.shape[2]):
1241     for i3 in range(a.shape[3]):
1242     z=""
1243     for i99 in range(a.shape[4]):
1244     if not z=="": z+="+"
1245     if o=="1":
1246 gross 439 z+="(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99]+b[i0,i1,i2,i3,i99],x,i99)
1247 gross 437 else:
1248 gross 439 z+="(%s)*%s[%s]**o+(%s)*%s[%s]"%(a[i0,i1,i2,i3,i99],x,i99,b[i0,i1,i2,i3,i99],x,i99)
1249 gross 437
1250     out+=tap+"%s[%s,%s,%s,%s]=%s\n"%(arg,i0,i1,i2,i3,z)
1251     return out
1252    
1253     def unrollLoopsOfGrad(a,b,o,arg,tap=""):
1254     out=""
1255     if a.rank==1:
1256     for i99 in range(a.shape[0]):
1257     if o=="1":
1258     out+=tap+"%s[%s]=(%s)\n"%(arg,i99,a[i99]+b[i99])
1259     else:
1260     out+=tap+"%s[%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i99,a[i99],i99,b[i99])
1261    
1262     elif a.rank==2:
1263     for i0 in range(a.shape[0]):
1264     for i99 in range(a.shape[1]):
1265     if o=="1":
1266     out+=tap+"%s[%s,%s]=(%s)\n"%(arg,i0,i99,a[i0,i99]+b[i0,i99])
1267     else:
1268     out+=tap+"%s[%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i99,a[i0,i99],i99,b[i0,i99])
1269    
1270     elif a.rank==3:
1271     for i0 in range(a.shape[0]):
1272     for i1 in range(a.shape[1]):
1273     for i99 in range(a.shape[2]):
1274     if o=="1":
1275     out+=tap+"%s[%s,%s,%s]=(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99]+b[i0,i1,i99])
1276     else:
1277     out+=tap+"%s[%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i99,a[i0,i1,i99],i99,b[i0,i1,i99])
1278    
1279     elif a.rank==4:
1280     for i0 in range(a.shape[0]):
1281     for i1 in range(a.shape[1]):
1282     for i2 in range(a.shape[2]):
1283     for i99 in range(a.shape[3]):
1284     if o=="1":
1285     out+=tap+"%s[%s,%s,%s,%s]=(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99]+b[i0,i1,i2,i99])
1286     else:
1287     out+=tap+"%s[%s,%s,%s,%s]=o*(%s)*x_ref[%s]**(o-1)+(%s)\n"%(arg,i0,i1,i2,i99,a[i0,i1,i2,i99],i99,b[i0,i1,i2,i99])
1288 gross 438 return out
1289 gross 441 def unrollLoopsOfDiv(a,b,o,arg,tap=""):
1290     out=tap+arg+"="
1291     if o=="1":
1292     z=0.
1293     for i99 in range(a.shape[0]):
1294     z+=b[i99,i99]+a[i99,i99]
1295     out+="(%s)"%z
1296     else:
1297     z=0.
1298     for i99 in range(a.shape[0]):
1299     z+=b[i99,i99]
1300     if i99>0: out+="+"
1301     out+="o*(%s)*x_ref[%s]**(o-1)"%(a[i99,i99],i99)
1302     out+="+(%s)"%z
1303     return out
1304    
1305 gross 437 def unrollLoopsOfInteriorIntegral(a,b,where,arg,tap=""):
1306     if where=="Function":
1307     xfac_o=1.
1308     xfac_op=0.
1309     z_fac=1./2.
1310     z_fac_s=""
1311     zo_fac_s="/(o+1.)"
1312     elif where=="FunctionOnBoundary":
1313     xfac_o=1.
1314     xfac_op=0.
1315     z_fac=1.
1316     z_fac_s="*dim"
1317     zo_fac_s="*(1+2.*(dim-1.)/(o+1.))"
1318     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1319     xfac_o=0.
1320     xfac_op=1.
1321     z_fac=1./2.
1322     z_fac_s=""
1323     zo_fac_s="/(o+1.)"
1324     out=""
1325     if a.rank==1:
1326     zo=0.
1327     zop=0.
1328     z=0.
1329     for i99 in range(a.shape[0]):
1330     if i99==0:
1331     zo+= xfac_o*a[i99]
1332     zop+= xfac_op*a[i99]
1333     else:
1334     zo+=a[i99]
1335     z+=b[i99]
1336    
1337     out+=tap+"%s=(%s)%s+(%s)%s"%(arg,zo,zo_fac_s,z*z_fac,z_fac_s)
1338     if zop==0.:
1339     out+="\n"
1340     else:
1341     out+="+(%s)*0.5**o\n"%zop
1342     elif a.rank==2:
1343     for i0 in range(a.shape[0]):
1344     zo=0.
1345     zop=0.
1346     z=0.
1347     for i99 in range(a.shape[1]):
1348     if i99==0:
1349     zo+= xfac_o*a[i0,i99]
1350     zop+= xfac_op*a[i0,i99]
1351     else:
1352     zo+=a[i0,i99]
1353     z+=b[i0,i99]
1354    
1355     out+=tap+"%s[%s]=(%s)%s+(%s)%s"%(arg,i0,zo,zo_fac_s,z*z_fac,z_fac_s)
1356     if zop==0.:
1357     out+="\n"
1358     else:
1359     out+="+(%s)*0.5**o\n"%zop
1360     elif a.rank==3:
1361     for i0 in range(a.shape[0]):
1362     for i1 in range(a.shape[1]):
1363     zo=0.
1364     zop=0.
1365     z=0.
1366     for i99 in range(a.shape[2]):
1367     if i99==0:
1368     zo+= xfac_o*a[i0,i1,i99]
1369     zop+= xfac_op*a[i0,i1,i99]
1370     else:
1371     zo+=a[i0,i1,i99]
1372     z+=b[i0,i1,i99]
1373    
1374     out+=tap+"%s[%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,zo,zo_fac_s,z*z_fac,z_fac_s)
1375     if zop==0.:
1376     out+="\n"
1377     else:
1378     out+="+(%s)*0.5**o\n"%zop
1379     elif a.rank==4:
1380     for i0 in range(a.shape[0]):
1381     for i1 in range(a.shape[1]):
1382     for i2 in range(a.shape[2]):
1383     zo=0.
1384     zop=0.
1385     z=0.
1386     for i99 in range(a.shape[3]):
1387     if i99==0:
1388     zo+= xfac_o*a[i0,i1,i2,i99]
1389     zop+= xfac_op*a[i0,i1,i2,i99]
1390    
1391     else:
1392     zo+=a[i0,i1,i2,i99]
1393     z+=b[i0,i1,i2,i99]
1394    
1395     out+=tap+"%s[%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,zo,zo_fac_s,z*z_fac,z_fac_s)
1396     if zop==0.:
1397     out+="\n"
1398     else:
1399     out+="+(%s)*0.5**o\n"%zop
1400    
1401     elif a.rank==5:
1402     for i0 in range(a.shape[0]):
1403     for i1 in range(a.shape[1]):
1404     for i2 in range(a.shape[2]):
1405     for i3 in range(a.shape[3]):
1406     zo=0.
1407     zop=0.
1408     z=0.
1409     for i99 in range(a.shape[4]):
1410     if i99==0:
1411     zo+= xfac_o*a[i0,i1,i2,i3,i99]
1412     zop+= xfac_op*a[i0,i1,i2,i3,i99]
1413    
1414     else:
1415     zo+=a[i0,i1,i2,i3,i99]
1416     z+=b[i0,i1,i2,i3,i99]
1417     out+=tap+"%s[%s,%s,%s,%s]=(%s)%s+(%s)%s"%(arg,i0,i1,i2,i3,zo,zo_fac_s,z*z_fac,z_fac_s)
1418     if zop==0.:
1419     out+="\n"
1420     else:
1421     out+="+(%s)*0.5**o\n"%zop
1422    
1423     return out
1424 gross 443 def unrollLoopsSimplified(b,arg,tap=""):
1425     out=""
1426     if isinstance(b,float) or b.rank==0:
1427     out+=tap+"%s=(%s)*x[0]\n"%(arg,str(b))
1428    
1429     elif b.rank==1:
1430     for i0 in range(b.shape[0]):
1431     out+=tap+"%s[%s]=(%s)*x[%s]\n"%(arg,i0,b[i0],i0)
1432     elif b.rank==2:
1433     for i0 in range(b.shape[0]):
1434     for i1 in range(b.shape[1]):
1435     out+=tap+"%s[%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,b[i0,i1],i1)
1436     elif b.rank==3:
1437     for i0 in range(b.shape[0]):
1438     for i1 in range(b.shape[1]):
1439     for i2 in range(b.shape[2]):
1440     out+=tap+"%s[%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,b[i0,i1,i2],i2)
1441     elif b.rank==4:
1442     for i0 in range(b.shape[0]):
1443     for i1 in range(b.shape[1]):
1444     for i2 in range(b.shape[2]):
1445     for i3 in range(b.shape[3]):
1446     out+=tap+"%s[%s,%s,%s,%s]=(%s)*x[%s]\n"%(arg,i0,i1,i2,i3,b[i0,i1,i2,i3],i3)
1447     return out
1448    
1449     def unrollLoopsOfL2(b,where,arg,tap=""):
1450     out=""
1451     z=[]
1452     if isinstance(b,float) or b.rank==0:
1453     z.append(b**2)
1454     elif b.rank==1:
1455     for i0 in range(b.shape[0]):
1456     z.append(b[i0]**2)
1457     elif b.rank==2:
1458     for i1 in range(b.shape[1]):
1459     s=0
1460     for i0 in range(b.shape[0]):
1461     s+=b[i0,i1]**2
1462     z.append(s)
1463     elif b.rank==3:
1464     for i2 in range(b.shape[2]):
1465     s=0
1466     for i0 in range(b.shape[0]):
1467     for i1 in range(b.shape[1]):
1468     s+=b[i0,i1,i2]**2
1469     z.append(s)
1470    
1471     elif b.rank==4:
1472     for i3 in range(b.shape[3]):
1473     s=0
1474     for i0 in range(b.shape[0]):
1475     for i1 in range(b.shape[1]):
1476     for i2 in range(b.shape[2]):
1477     s+=b[i0,i1,i2,i3]**2
1478     z.append(s)
1479     if where=="Function":
1480     xfac_o=1.
1481     xfac_op=0.
1482     z_fac_s=""
1483     zo_fac_s=""
1484     zo_fac=1./3.
1485     elif where=="FunctionOnBoundary":
1486     xfac_o=1.
1487     xfac_op=0.
1488     z_fac_s="*dim"
1489     zo_fac_s="*(2.*dim+1.)/3."
1490     zo_fac=1.
1491     elif where in ["FunctionOnContactZero","FunctionOnContactOne"]:
1492     xfac_o=0.
1493     xfac_op=1.
1494     z_fac_s=""
1495     zo_fac_s=""
1496     zo_fac=1./3.
1497     zo=0.
1498     zop=0.
1499     for i99 in range(len(z)):
1500     if i99==0:
1501     zo+=xfac_o*z[i99]
1502     zop+=xfac_op*z[i99]
1503     else:
1504     zo+=z[i99]
1505     out+=tap+"%s=sqrt((%s)%s"%(arg,zo*zo_fac,zo_fac_s)
1506     if zop==0.:
1507     out+=")\n"
1508     else:
1509     out+="+(%s))\n"%(zop*0.5**2)
1510     return out
1511 gross 493 #=======================================================================================================
1512     # transpose
1513     #=======================================================================================================
1514     def transposeTest(r,offset):
1515     if isinstance(r,float): return r
1516     s=r.shape
1517     s1=1
1518     for i in s[:offset]: s1*=i
1519     s2=1
1520     for i in s[offset:]: s2*=i
1521     out=numarray.reshape(r,(s1,s2))
1522     out.transpose()
1523     return numarray.resize(out,s[offset:]+s[:offset])
1524 gross 443
1525 gross 493 name,tt="transpose",transposeTest
1526     for case0 in ["array","Symbol","constData","taggedData","expandedData"]:
1527     for sh0 in [ (), (3,), (4,5), (6,2,2),(3,2,3,4)]:
1528     for offset in range(len(sh0)+1):
1529     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1530     tname="test_%s_%s_rank%s_offset%s"%(name,case0,len(sh0),offset)
1531     text+=" def %s(self):\n"%tname
1532     sh_t=sh0[offset:]+sh0[:offset]
1533    
1534     # sh_t=list(sh0)
1535     # sh_t[offset+1]=sh_t[offset]
1536     # sh_t=tuple(sh_t)
1537     # sh_r=[]
1538     # for i in range(offset): sh_r.append(sh0[i])
1539     # for i in range(offset+2,len(sh0)): sh_r.append(sh0[i])
1540     # sh_r=tuple(sh_r)
1541    
1542     a_0=makeArray(sh0,[-1.,1])
1543     if case0 in ["taggedData", "expandedData"]:
1544     a1_0=makeArray(sh0,[-1.,1])
1545     else:
1546     a1_0=a_0
1547     r=tt(a_0,offset)
1548     r1=tt(a1_0,offset)
1549     text+=mkText(case0,"arg",a_0,a1_0)
1550     text+=" res=%s(arg,%s)\n"%(name,offset)
1551     if case0=="Symbol":
1552     text+=mkText("array","s",a_0,a1_0)
1553     text+=" sub=res.substitute({arg:s})\n"
1554     res="sub"
1555     text+=mkText("array","ref",r,r1)
1556     else:
1557     res="res"
1558     text+=mkText(case0,"ref",r,r1)
1559     text+=mkTypeAndShapeTest(case0,sh_t,"res")
1560     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1561    
1562     if case0 == "taggedData" :
1563     t_prog_with_tags+=text
1564     else:
1565     t_prog+=text
1566    
1567     print test_header
1568     # print t_prog
1569     print t_prog_with_tags
1570     print test_tail
1571     1/0
1572 gross 441 #=======================================================================================================
1573 gross 443 # L2
1574     #=======================================================================================================
1575     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1576     for data in ["Data","Symbol"]:
1577     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1578     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1579     tname="test_L2_on%s_from%s_rank%s"%(where,data,len(sh))
1580     text+=" def %s(self):\n"%tname
1581     text+=" \"\"\"\n"
1582     text+=" tests L2-norm of %s on the %s\n\n"%(data,where)
1583     text+=" assumptions: self.domain supports integration on %s\n"%where
1584     text+=" \"\"\"\n"
1585     text+=" dim=self.domain.getDim()\n"
1586     text+=" w=%s(self.domain)\n"%where
1587     text+=" x=w.getX()\n"
1588     o="1"
1589     if len(sh)>0:
1590     sh_2=sh[:len(sh)-1]+(2,)
1591     sh_3=sh[:len(sh)-1]+(3,)
1592     b_2=makeArray(sh[:len(sh)-1]+(2,),[-1.,1])
1593     b_3=makeArray(sh[:len(sh)-1]+(3,),[-1.,1])
1594     else:
1595     sh_2=()
1596     sh_3=()
1597     b_2=makeArray(sh,[-1.,1])
1598     b_3=makeArray(sh,[-1.,1])
1599    
1600     if data=="Symbol":
1601     val="s"
1602     res="sub"
1603     else:
1604     val="arg"
1605     res="res"
1606     text+=" if dim==2:\n"
1607     if data=="Symbol":
1608     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_2)
1609    
1610     text+=" %s=Data(0,%s,w)\n"%(val,sh_2)
1611     text+=unrollLoopsSimplified(b_2,val,tap=" ")
1612     text+=unrollLoopsOfL2(b_2,where,"ref",tap=" ")
1613     text+="\n else:\n"
1614     if data=="Symbol":
1615     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh_3)
1616     text+=" %s=Data(0,%s,w)\n"%(val,sh_3)
1617     text+=unrollLoopsSimplified(b_3,val,tap=" ")
1618     text+=unrollLoopsOfL2(b_3,where,"ref",tap=" ")
1619     text+="\n res=L2(arg)\n"
1620     if data=="Symbol":
1621     text+=" sub=res.substitute({arg:s})\n"
1622     text+=" self.failUnless(isinstance(res,Symbol),\"wrong type of result.\")\n"
1623     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1624     else:
1625     text+=" self.failUnless(isinstance(res,float),\"wrong type of result.\")\n"
1626     text+=" self.failUnlessAlmostEqual(%s,ref,int(-log10(self.RES_TOL)),\"wrong result\")\n"%res
1627     t_prog+=text
1628     print t_prog
1629     1/0
1630    
1631     #=======================================================================================================
1632 gross 441 # div
1633     #=======================================================================================================
1634     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1635     for data in ["Data","Symbol"]:
1636     for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1637     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1638     tname="test_div_on%s_from%s_%s"%(where,data,case)
1639     text+=" def %s(self):\n"%tname
1640     text+=" \"\"\"\n"
1641     text+=" tests divergence of %s on the %s\n\n"%(data,where)
1642     text+=" assumptions: %s(self.domain) exists\n"%case
1643     text+=" self.domain supports div on %s\n"%where
1644     text+=" \"\"\"\n"
1645     if case=="ReducedSolution":
1646     text+=" o=1\n"
1647     o="1"
1648     else:
1649     text+=" o=self.order\n"
1650     o="o"
1651     text+=" dim=self.domain.getDim()\n"
1652     text+=" w_ref=%s(self.domain)\n"%where
1653     text+=" x_ref=w_ref.getX()\n"
1654     text+=" w=%s(self.domain)\n"%case
1655     text+=" x=w.getX()\n"
1656     a_2=makeArray((2,2),[-1.,1])
1657     b_2=makeArray((2,2),[-1.,1])
1658     a_3=makeArray((3,3),[-1.,1])
1659     b_3=makeArray((3,3),[-1.,1])
1660     if data=="Symbol":
1661     text+=" arg=Symbol(shape=(dim,),dim=dim)\n"
1662     val="s"
1663     res="sub"
1664     else:
1665     val="arg"
1666     res="res"
1667     text+=" %s=Vector(0,w)\n"%val
1668     text+=" if dim==2:\n"
1669     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1670     text+=unrollLoopsOfDiv(a_2,b_2,o,"ref",tap=" ")
1671     text+="\n else:\n"
1672    
1673     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1674     text+=unrollLoopsOfDiv(a_3,b_3,o,"ref",tap=" ")
1675     text+="\n res=div(arg,where=w_ref)\n"
1676     if data=="Symbol":
1677     text+=" sub=res.substitute({arg:s})\n"
1678     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1679     text+=" self.failUnlessEqual(res.getShape(),(),\"wrong shape of result.\")\n"
1680     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong function space of result.\")\n"%res
1681     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1682 gross 437
1683 gross 441
1684     t_prog+=text
1685     print t_prog
1686     1/0
1687    
1688 gross 439 #=======================================================================================================
1689     # interpolation
1690     #=======================================================================================================
1691     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne","Solution","ReducedSolution"]:
1692     for data in ["Data","Symbol"]:
1693     for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1694     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1695     if where==case or \
1696     ( case in ["ContinuousFunction","Solution","ReducedSolution"] and where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"] ) or \
1697     ( case in ["FunctionOnContactZero","FunctionOnContactOne"] and where in ["FunctionOnContactZero","FunctionOnContactOne"] ) or \
1698     (case=="ContinuousFunction" and where in ["Solution","ReducedSolution"]) or \
1699     (case=="Solution" and where=="ReducedSolution") :
1700    
1701    
1702     text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1703     tname="test_interpolation_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1704     text+=" def %s(self):\n"%tname
1705     text+=" \"\"\"\n"
1706     text+=" tests interpolation for rank %s %s onto the %s\n\n"%(len(sh),data,where)
1707     text+=" assumptions: self.domain supports inpterpolation from %s onto %s\n"%(case,where)
1708     text+=" \"\"\"\n"
1709     if case=="ReducedSolution" or where=="ReducedSolution":
1710     text+=" o=1\n"
1711     o="1"
1712     else:
1713     text+=" o=self.order\n"
1714     o="o"
1715     text+=" dim=self.domain.getDim()\n"
1716     text+=" w_ref=%s(self.domain)\n"%where
1717     text+=" x_ref=w_ref.getX()\n"
1718     text+=" w=%s(self.domain)\n"%case
1719     text+=" x=w.getX()\n"
1720     a_2=makeArray(sh+(2,),[-1.,1])
1721     b_2=makeArray(sh+(2,),[-1.,1])
1722     a_3=makeArray(sh+(3,),[-1.,1])
1723     b_3=makeArray(sh+(3,),[-1.,1])
1724     if data=="Symbol":
1725     text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1726     val="s"
1727     res="sub"
1728     else:
1729     val="arg"
1730     res="res"
1731     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1732     text+=" ref=Data(0,%s,w_ref)\n"%str(sh)
1733     text+=" if dim==2:\n"
1734     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1735     text+=unrollLoops(a_2,b_2,o,"ref",tap=" ",x="x_ref")
1736     text+=" else:\n"
1737    
1738     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1739     text+=unrollLoops(a_3,b_3,o,"ref",tap=" ",x="x_ref")
1740     text+=" res=interpolate(arg,where=w_ref)\n"
1741     if data=="Symbol":
1742     text+=" sub=res.substitute({arg:s})\n"
1743     text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1744     text+=" self.failUnlessEqual(%s.getFunctionSpace(),w_ref,\"wrong functionspace of result.\")\n"%res
1745     text+=" self.failUnlessEqual(res.getShape(),%s,\"wrong shape of result.\")\n"%str(sh)
1746     text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1747     t_prog+=text
1748     print test_header
1749     print t_prog
1750     print test_tail
1751     1/0
1752 gross 437
1753 gross 429 #=======================================================================================================
1754 gross 438 # grad
1755 gross 437 #=======================================================================================================
1756     for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1757     for data in ["Data","Symbol"]:
1758 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution"]:
1759     for sh in [ (),(2,), (4,5), (6,2,2)]:
1760 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1761 gross 438 tname="test_grad_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1762 gross 437 text+=" def %s(self):\n"%tname
1763     text+=" \"\"\"\n"
1764 gross 438 text+=" tests gradient for rank %s %s on the %s\n\n"%(len(sh),data,where)
1765 gross 437 text+=" assumptions: %s(self.domain) exists\n"%case
1766 gross 438 text+=" self.domain supports gardient on %s\n"%where
1767 gross 437 text+=" \"\"\"\n"
1768     if case=="ReducedSolution":
1769     text+=" o=1\n"
1770     o="1"
1771     else:
1772     text+=" o=self.order\n"
1773     o="o"
1774     text+=" dim=self.domain.getDim()\n"
1775     text+=" w_ref=%s(self.domain)\n"%where
1776 gross 438 text+=" x_ref=w_ref.getX()\n"
1777 gross 437 text+=" w=%s(self.domain)\n"%case
1778     text+=" x=w.getX()\n"
1779     a_2=makeArray(sh+(2,),[-1.,1])
1780     b_2=makeArray(sh+(2,),[-1.,1])
1781     a_3=makeArray(sh+(3,),[-1.,1])
1782     b_3=makeArray(sh+(3,),[-1.,1])
1783     if data=="Symbol":
1784 gross 438 text+=" arg=Symbol(shape=%s,dim=dim)\n"%str(sh)
1785 gross 437 val="s"
1786     res="sub"
1787     else:
1788     val="arg"
1789     res="res"
1790     text+=" %s=Data(0,%s,w)\n"%(val,str(sh))
1791 gross 438 text+=" ref=Data(0,%s+(dim,),w_ref)\n"%str(sh)
1792 gross 437 text+=" if dim==2:\n"
1793     text+=unrollLoops(a_2,b_2,o,val,tap=" ")
1794 gross 438 text+=unrollLoopsOfGrad(a_2,b_2,o,"ref",tap=" ")
1795 gross 437 text+=" else:\n"
1796    
1797     text+=unrollLoops(a_3,b_3,o,val,tap=" ")
1798 gross 438 text+=unrollLoopsOfGrad(a_3,b_3,o,"ref",tap=" ")
1799     text+=" res=grad(arg,where=w_ref)\n"
1800 gross 437 if data=="Symbol":
1801     text+=" sub=res.substitute({arg:s})\n"
1802 gross 438 text+=" self.failUnless(isinstance(res,%s),\"wrong type of result.\")\n"%data
1803     text+=" self.failUnlessEqual(res.getShape(),%s+(dim,),\"wrong shape of result.\")\n"%str(sh)
1804 gross 437 text+=" self.failUnless(Lsup(%s-ref)<=self.RES_TOL*Lsup(ref),\"wrong result\")\n"%res
1805    
1806    
1807     t_prog+=text
1808     print test_header
1809     print t_prog
1810     print test_tail
1811     1/0
1812 gross 438
1813    
1814 gross 437 #=======================================================================================================
1815 gross 438 # integrate
1816 gross 437 #=======================================================================================================
1817 gross 438 for where in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1818 gross 437 for data in ["Data","Symbol"]:
1819 gross 438 for case in ["ContinuousFunction","Solution","ReducedSolution","Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]:
1820     for sh in [ (),(2,), (4,5), (6,2,2),(4,5,3,2)]:
1821     if (not case in ["Function","FunctionOnBoundary","FunctionOnContactZero","FunctionOnContactOne"]) or where==case:
1822 gross 437 text=" #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
1823 gross 438 tname="test_integrate_on%s_from%s_%s_rank%s"%(where,data,case,len(sh))
1824 gross 437 text+=" def %s(self):\n"%tname
1825     text+=" \"\"\"\n"
1826 gross 438 text+=" tests integral of rank %s %s on the %s\n\n"%(len(sh),data,where)
1827     text+=" assumptions: %s(self.domain) exists\n"%case
1828     text+=" self.domain supports integral on %s\n"%where
1829    
1830 gross 437 text+=" \"\"\"\n"
1831     if case=="ReducedSolution":
1832     text+=" o=1\n"
1833     o="1"
1834     else: