Telecomms 3D Erlang Blocking Problem-Efficiency improvements please

*To*: mathgroup at smc.vnet.net*Subject*: [mg19512] Telecomms 3D Erlang Blocking Problem-Efficiency improvements please*From*: king at dircon.co.uk (Nigel King)*Date*: Sat, 28 Aug 1999 15:53:28 -0400*Sender*: owner-wri-mathgroup at wolfram.com

Dear Mathematica Group, I am working on a Telecoms Erlang blocking problem which uses a large 3 dimensional probability array which needs evaluating for 58x58x38 if possible. Unfortunately I have only got to 18x18x12 so far before running out of memory (64Mbyte Kernel) and the time is rising very high. I have done one of these problmes before which was 6 dimensions but there was symmetry which could be exploited enabling a C program to be written for the central summations. In this case there are no symmetries (that I can see) so the probability array is 'Solved' using simultaneous equations. The core program below seems to me to be elegant, and I validated that it produces the correct results for simple cases and that the complex cases are within a few percent of the result which I would expect. However, I need if possible to improve the capability of the program. I shall obviously need to experiment with better methods but would appreciate some pointers. Here are some thoughts which I have had:- 1. P[n1,n2,n3] is the probability array, should I have used P[[n1,n2,n3]] 2. Solve should be replaced by NSolve. I have tried this but the pointers n1, n2 and n3 must be integers. Is there a way to force N[] not to wreck the array pointers. I wondered about evaluating ToExpession["P"<>ToString[n1]<>"n"<>ToString[n2]<>"n"<>ToString[n3]] for the variables. 3. Re-reading "Power Programing with Mathematica The Kernel" 4. Any ideas from this group. Thanks in advance Nigel King The program follows P[] is the probability matrix Peq[] are the balance equations \[Lambda] are birth flows which have many conditions \[Mu] are deaths x is the set of equations with one dropped because of over constraint y is all the terms of the probability matrix tot is the sum of all probabilities (for normalising the result) b1, b2, bt and bs are the blocking probabilities Parray1 is the Program which starts by clearing the previous results c1, c2 and c3 are numbers of circuits a result with timing In[143]:= Timing[Parray[.1, .1, .1, .1, 3, 3, 3]] Out[143]= \!\({2.0500000000029104`\ Second, {1.2327646535342597`*^-7, 1.2327646535342602`*^-7, 0.0001512386148015114`, 2.9825151264562333`*^-10, 0.00003787136649593273`}}\) Core program below had to change "less than or equal symbol" to <= so may not paste into Mathematica without reconverting MacSoup complained about not Latin Characterset Parray1[a1_, a2_, a3_, a4_, d1_, d2_, d3_] := (Clear[Peq, P, \[Lambda], as, \[Mu], a, c1, c2, c3]; c1 = d1; c2 = d2; c3 = d3; Peq[n1_, n2_, n3_] := P[n1, n2, n3] == (\[Mu]((n1 + 1)P[n1 + 1, n2, n3] + (n2 + 1)P[n1, n2 + 1, n3] + (n3 + 1)P[n1, n2, n3 + 1]) + P[n1 - 1, n2, n3]\[Lambda][1, n1 - 1, n2, n3] + P[n1, n2 - 1, n3]\[Lambda][2, n1, n2 - 1, n3] + P[n1, n2, n3 - 1]\[Lambda][3, n1, n2, n3 - 1])/(\[Lambda][1, n1, n2, n3] + \[Lambda][2, n1, n2, n3] + \[Lambda][3, n1, n2, n3] + \[Mu](n1 + n2 + n3)); P[0, 0, 0] = 1; P[n1_, n2_, n3_] /; n1 > c1 || n2 > c2 || n3 > c3 || n1 < 0 || n2 < 0 || n3 < 0 := P1[n1, n2, n3] = 0; \[Mu] = 1; \[Lambda][1, n1_, n2_, n3_] := If[(0 <= n1 < c1) && (0 <= n2 <= c2) && (0 <= n3 <= c3), a1 + a4(c1 - n1)/(c1 + c2 - n1 - n2), 0]; \[Lambda][2, n1_, n2_, n3_] := If[(0 <= n1 <= c1) && (0 <= n2 < c2) && (0 <= n3 <= c3), a2 + a4(c2 - n2)/(c1 + c2 - n1 - n2), 0]; \[Lambda][3, n1_, n2_, n3_] := If[(n1 == c1) && (n2 == c2) && (0 <= n3 < c3), a1 + a2 + a3 + a4, If[(n1 == c1) && (0 <= n2 < c2) && (0 <= n3 < c3), a1 + a3, If[(0 <= n1 < c1) && (n2 == c2) && (0 <= n3 < c3), a2 + a3, If[(0 <= n1 < c1) && (0 <= n2 < c2) && (0 <= n3 < c3), a3, 0]]]]; x = Simplify[ Drop[Flatten[ Table[Peq[n1, n2, n3], {n1, 0, c1}, {n2, 0, c2}, {n3, 0, c3}]], 1]]; y = Drop[Complement[ Flatten[y1 = Table[P[n1, n2, n3], {n1, 0, c1}, {n2, 0, c2}, {n3, 0, c3}]], {1}], 0]; a = Simplify[Flatten[Solve[x, y, y]]]; tot = Apply[Plus, y] + 1; b1 = Apply[Plus, Flatten[Table[P[c1, n2, c3], {n2, 0, c2}]]]/tot; b2 = Apply[Plus, Flatten[Table[P[n1, c2, c3], {n1, 0, c1}]]]/tot; bt = Apply[Plus, Flatten[Table[P[n1, n2, c3], {n1, 0, c1}, {n2, 0, c2}]]]/ tot; bs = P[c1, c2, c3]/tot;);