Carl Love

Carl Love

28070 Reputation

25 Badges

13 years, 31 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@acer wrote

  • Maple doesn't offer precomputing the P,LU for a band storage, in part because it needs a structure larger than input AT and things get a bit awkward functionality-wise.

I just wrote a triadiagonal solver that stores prefactors and thus beats that NAG band solver from your example by a factor of 5 timewise. The prefactors could be stored in place in the tridiagonal Matrix (without altering the fact that it is tridiagonal), but I didn't implement that aspect. Indeed, they could simply replace the main diagonal and 1st subdiagonal.

In my example, rather than solving repeatedly using the same side vector, I use all different side vectors.
 

restart:

N:= 2^12: #space dimension size
maxiter:= 2^7: #time dimension size

AT:= LinearAlgebra:-RandomMatrix(
   N$2, generator= 0.0..100.0, datatype= hfloat, shape= band[1,1]
):

V:= LinearAlgebra:-RandomMatrix((N,maxiter), generator= -100.0..100.0, datatype= hfloat):
XT:= Matrix((N,maxiter), datatype= hfloat):

#Time and accuracy test of Maple's/NAG's TriD solver:

st:= time():
for i to maxiter do
   XT[..,i]:= LinearAlgebra:-LinearSolve(AT,V[..,i])
end do:
time()-st;
LinearAlgebra:-Norm(AT.XT - V);

.172

0.210472710371334415e-7

#Carl's 2-stage TriD solver:
TriD_prefactor_inner:= proc(
   A::Array(datatype= float[8]),
   B::Array(datatype= float[8]),
   C::Array(datatype= float[8]),
   n::integer[4]
)::integer[4];
local k::integer[4];
   for k to n do A[k]:= A[k]/B[k]; B[k+1]:= B[k+1] - A[k]*C[k] od;
   0 #Just to return something
end proc:

TriD_prefactor_comp:= Compiler:-Compile(TriD_prefactor_inner):

TriD_prefactor:= proc(M::Matrix(square, storage= band[1,1], datatype= hfloat))
local   
   n:= rhs(rtable_dims(M)[1])-1,
   A:= rtable(1..n, k-> M[k+1, k], datatype= hfloat),
   B:= rtable(1..n+1, k-> M[k, k], datatype= hfloat),
   C:= rtable(1..n, k-> M[k, k+1], datatype= hfloat)
;
   TriD_prefactor_comp(A, B, C, n);
   (A,B,C)
end proc:

TriD_solve_inner:= proc(
   A::Array(datatype= float[8]),
   B::Array(datatype= float[8]),
   C::Array(datatype= float[8]),
   X::Array(datatype= float[8]),
   n::integer[4]
)::integer[4];
local k::integer[4];
   for k to n do X[k+1]:= X[k+1] - A[k]*X[k] od;
   X[n+1]:= X[n+1]/B[n+1];
   for k from n by -1 to 1 do X[k]:= (X[k]-C[k]*X[k+1])/B[k] od;
   0 #Just to return something
end proc:

TriD_solve_comp:= Compiler:-Compile(TriD_solve_inner):

TriD_solve:= proc(A, B, C, D):
local n:= rtable_num_elems(A), X:= rtable(D);
   TriD_solve_comp(A, B, C, X, n);
   X
end proc:
   

XT:= Matrix((N,maxiter), datatype= hfloat): #Clear old solutions.

#Time and accuracy test of Carl's 2-stage solver
st:= time():
AT_pf:= TriD_prefactor(AT):
for i to maxiter do
   XT[..,i]:= TriD_solve(AT_pf, V[..,i])
end do:
time()-st;
LinearAlgebra:-Norm(AT.XT - V);

0.32e-1

0.468513432494432891e-7

 


 

Download TriD_prefactored.mw

@torabi  GIGO[*1]: Translations can only work if the original is correct. My Matlab knowledge is quite rusty, but in this case, it looks like must be a Matrix for your code to make any sense in either language. Yet, you have M:= 2, making a scalar. So, I must ask, Did you ever run the code in Matlab? And did it work?

[*1]GIGO stands for "Garbage In, Garbage Out", a tongue-in-cheek law of computer science.

Vote up for the derivation and exposition.

Regarding the iterative solver: For n x tridiagonal A, the system A.U = B can be solved in O(n) steps and O(1) auxilliary memory, which is much faster than even multiplying by A^(-1), let alone computing A^(-1) in the first place. And often the memory alone needed for A^(-1) is prohibitive. See the Wikipedia article "Tridiagonal matrix algorithm".

It works fine, like this:

IfProc:= x-> if x=0 then 1 elif x=1 then 2 elif x=3 then 7 else 10 fi:
IfProc~([$0..4]);
                       [1, 2, 10, 7, 10]

 

@shayas75 I'm not sure what you mean. If you type exactly what I entered in red (the color makes no difference), you should get the answer 2. Isn't that good enough?

@torabi 

Matlab           Maple
-----------------------------------------------------
Q(.., k)         Q[.., k]  
size(M,k)        (1+rhs-lhs)([rtable_dims(M)][k])
T= S:h:F         T:= <seq(S..F, h)>
length(T)        1+max((rhs-lhs)~([rtable_dims(T)]))
zeros(N,lt)      Matrix(N,lt)
A\D              LinearAlgebra:-LinearSolve(A,D) 
for n = M:N      for n from N to M do
for n = M:S:N    for n from N by S to M do 
for n = [a,b,c]  for n in [a,b,c] do

Anticipating your next question, I also translated some for loop syntax.

Here's a detailed example of finding the conjugate of a partition via its Ferrers diagram. The diagram technique verifies the correspondence between partitions and their conjugates and also provides an effective technique for computing the conjugate.

Consider the partition: 14 = 1 + 3 + 4 + 6. It has length 4 and maximum 6. The conjugate will have length 6 and maximum 4. Here's the Ferrers diagram:

6: @ @ @ @ @ @
4: @ @ @ @ 
3: @ @ @
1: @
----------------------------Now count the entries in each column
   4 3 3 2 1 1

So, the conjugate partition is 14 = 1 + 1 + 2 + 3 + 3 + 4. Computers are better at adding 1s than counting @s, so in my program, it's like this:

  [0, 0, 0, 0, 0, 1] +  #1
  [0, 0, 0, 1, 1, 1] +  #3
  [0, 0, 1, 1, 1, 1] +  #4
  [1, 1, 1, 1, 1, 1]    #6 
----------------------------Add these as vectors
  [1, 1, 2, 3, 3, 4]

Now suppose that you called my program with arguments (20, 6). It would first call combinat:-partition(14, 6). One of the partitions in the returned list would be [1, 3, 4, 6]. The Ferrers process in my program converts that to [1, 1, 2, 3, 3, 4]. Finally, 1 is added to each entry, yielding [2, 2, 3, 4, 4, 5], which is a partition of 20 of length 6, which is what was desired.

@ms0439883 I said that the min(n-k, k) was needed for the case k > n/2. In your example, n=5, k=2, so k < n/2. The min is only required to satisfy the overly zealous argument checker in combinat:-partition. It has nothing to do with the logic of the my algorithm, with conjugates of partitions, or with Ferrers diagrams. 

If you haven't yet, read that Wikipedia article that I linked, which'll explain how Ferrers diagrams are used to get the conjugate. Later, I'll post an example. 

If this Question is pursuant to your other Question about a system of 4 degree-3 polynomials with 4 variables and 3 symbolic parameters: The Newton-Raphson method cannot work with symbolic parameters.

I don't think that there's any way around it, in Maple or any other system: If you want usable solutions, then you'll need to supply numeric values for the parameters. If you do supply numeric values, the command fsolve should provide results very quickly, such that you could, for example, make plots of results versus parameters.

@brian bovril But you yourself said 

  • Just to clarify, each contestant was given an opaque bag containing a buff. Once they were dispensed, they were told to reveal their buffs (simultaneously). So that would be Scenario 2.

And aren't we all in agreement that the probability for Scenario 2 is 6/1001, whereas 1/1001 is the probability for Scenario 1?? (MMcdara used the word "strategy" instead of "scenario".)

@vv Oh my, I see that most of the new syntax features added for Maple 2019 are likewise not supported in 2-D input. I foresee this being a nuisance for those of us who answer questions here.

@ms0439883 

For easy reference, here is the line of code that I'm explaining:

[seq(1 +~ `+`(seq([0$(k-p), 1$p], p= P)), P= combinat:-partition(n-k, min(n-k, k)))]

Maple's combinat:-partition(n,k) returns a list of the partitions of n whose maximal elements are at most k. I need to conjugate those partitions whose maximal elements are equal to k, so I start with partitions of n-k whose maximal elements are at most k with the plan to compensate for exactly one missing k element later. The syntax of partition(n,k) won't allow k > n, so partition(n-k, k) would give an error for k > n/2, so I use partition(n-k, min(n-k, k)).

Now I have a list of partitions to process. The command [seq(f(x), x= L)] iterates over a list L, assigning each element to x in succession, passing it to f, and forming a list of the results. Hence my outer seq command is [seq(C(P), P= ...)] where is some as-yet-unexplained operation that will both conjugate P and compensate for the missing k. You must now remember that uppercase is always an individual partition of n-k expressed as a list of integers; P is neither a list of partitions nor an individual integer in a partition.

To conjugate a partition, I essentially use a Ferrers diagram. I'll assume that you either know what that is, or you can look it up (see, for example, this Wikipedia article). Instead of dots, I use 1s. For syntax reasons, I need lists of the same length, so I use lists of 1s padded on the left with 0s with each list having length k. The elementwise sum of these lists is the conjugate.

The Maple command X$N, where is a nonnegative integer, produces a sequence of copies of X. Thus [0$(k-p), 1$p] produces a length-list of p 1s padded on the left with 0s, exactly as described in the previous paragraph. Thus, seq([0$(k-p), 1$p], p= P) produces (as a sequence of lists) the entire Ferrers diagram of P. Note that lowercase p iterates over the integers in P.

Most of Maple's binary infix operators can be used in functional form (i.e., coming before their arguments) by enclosing them in backquotes. Thus `+`(5, 3, 7) is the same as 5+3+7. If the arguments are lists of the same length (as in the Ferrers diagram), then they are added elementwise. 

Most of Maple's binary infix operators can also be used in an infix elementwise form by appending tilde to the operator. Thus 1 +~ [x,y] is the same as [1+x, 1+y]. The purpose of this in my code is to compensate for the missing k by adding 1 to each element of each conjugated partition.

Please let me know if you understand all this or if you have further questions.

@Kalithiel The flag at the top of your Question says that you're using Maple 2019. The local syntax that I used will work in Maple 2019, but not in any earlier version. So, I think that you're using an earlier Maple. But, no matter, the old syntax (as shown by VV) will work in both new and old versions.

@nm If an item has Replies, then the Replies must be deleted first before the item itself can be deleted. By the way, you are a moderator.

I don't see why your Answer or any Answer should be deleted simply because it's not the most-complete or most-general answer possible. That would be a very high bar in many cases, and some things would get no Answers. Perhaps the OP does indeed only need the size-2 case.

@nm Of course, that only works when the size is 2. But I suspect that it could be applied recursively.

First 273 274 275 276 277 278 279 Last Page 275 of 709