@oldstudent There are several ways to ensure that the origin is included in the view.

For example, given your DataTable3,

plots:-display(plot([[0,0]]), plots:-pointplot( DataTable3 ));

minx,maxx:=(min,max)(DataTable3[..,1]):
plots:-pointplot( DataTable3, view=[min(-0.1,miny)..maxx, default] ); 

@oldstudent Do you want,

plots:-pointplot( DataTable3 );

That routine accepts an m-by-2 Matrix, and interprets the two columns as x and y values.

@jonlg 

• The help-page for topic rtable_indexing may be of use, if you haven't seen it yet. It ought to be linked from more pages than it is, IMO.
  
  (If I used round brackets to index into R, it was probably by omission, as I was editing your own second attachment. I try not to use round-bracket "programmer indexing" unless I need to, so that I'll know there was a reason for it when I return to the worksheet after a long time.)
  
  
• If you haven't read it yet, the Programming Manual is worth a read. (Download from here.)

@fereydoon_shekofte bytesused is not right, since it can report the bytes processed by garbage collection that may have nothing to do with the storage of the expression just assigned to a name. It is not a representation of the memory used to store anything in particular. The "used" in "bytesused" does not refer to memory that is storing a thing.

Up until recently the difference in kernelopts(bytesalloc) was a somewhat OK measure. But it can be way off, in modern versions, as garbage collection can now release memory to the OS. It was never perfect for this, either.

Part of the difficulty is that the operation to construct a thing can involve memory for temporary objects that are no longer present upon completion of the construction. It seems to me that you don't want that counted, is that right? You want only the bytes used to store a thing internally?

If you want a measure that is somewhat near proportional to the memory needed only to store a value then perhaps you could look at its `length`, or the `length` of it sprintf'd to %m format (so-called dotm). Those are not always in bytes, but they might be reasonably close to proportional. It could be tricky to handle modules.

As best I can guess at present, you might be wanting to do something like this...

MapleSimulation_test_a0.mw

@jonlg You've used `A` in two different ways. Let's keep the line that assigns to `A` the uninitialized 4-dimensional Array.

Let's focus on the first loop. You don't want to assign that operator (taking 4 parameters n,b,theta,m) to the name `A` as that will just clobber the earlier Array assignment. Let's assign that operator to the name `Afunc` instead, for now.

You may not need Afunc to be defined inside the loop. It'll depend on how this all gets unmuddled.

You can use square-bracket indexing for R, on the left hand side of that assignment in that first loop. It's makes things more clear (IMO) if you only use round-brackets for indexing when it is necessary. It's less confusing when you read this years later. I'll use you round-bracketing style below, even though I'm pretty such it mostly?! means indexing.

Now, a key thing: In your assignment to R(n+1) you have the term A(n). What is does that mean? As you had it as an operator `A` takes 4 parameters. I suspect that your A(n) here is intended as a function call to your operator `A`. But if so then what are the 2nd, 3rd, and 4th parameters (ie. b, theta, and m) supposed to be for that term A(n) in that assignment to R(n+1)?

Are you trying to set it up so that for each and every any fixed set of values of b, theta, and m the sequence of R values and A values should be computed using those formulae?

Could you be more precise than just, "generalized symbolic answer"? What kind of terms would an acceptable answer have? What do you know about omega[1] and omega[2], relative to each other?

acer

@H-R The command that does this is LinearAlgebra:-GenerateMatrix.

But it takes a fast machine about 3-4 minutesfor GenerateMatrix to process your system of 2500x2500 equations and variables. This is the wrong way to go. As I wrote before, much better would be to never form the symbolic equations at all, if you need to deal with problems of this size. You seem insistent on ignoring this key point. You should instead write code that writes numeric coefficients directly to the lhs Matrix M and the rhs Vector b, to represent M.x=b.

Your problem has more serious issues. The numeric coefficients range from 10^1000 down to -10^1000 and this will not fit inside the range of hardware double precision. Casting your data to float[8] data type will result in many entries becoming infinity, and that will not solve using LU, QR, or SVD. Moving from hardware to software working precision will make it run too long - expectedly, I could get no factorisation after even hours and hours... So, either your code which generates the system is just wrong, or your problem is very badly scaled, and your code doesn't correct t that. This is a serious problem with your code/approach.

Another issue is that the rhs of the M.x=b representation, b, seems to be almost all zero. So the system may be inconsistent, or under specified, or perhaps wrongly generated. Hard to test his serious this is until the scaling problem of entries of M is resolved.

Written well a FEM calculation producing 2500x2500/systems should be able to compute from start to end in 2-20 seconds or so on a modern machine, even within Maple. It seems that you have at least three serious problems which will either make that speed impossible or even make any solution unattainable.

@H-R In that case you should separate the unknowns from the coefficients, to produce a purely numeric Matrix M and a purely numeric Vector rhs b such that M.x=b would hold. Here x stands for the vector of unknowns in which your system is linear. For decent petformamce M must be floats not exact rationals, and preferaby have datatype=float[8].

I suspect that it would have been much better to never have formed your symbolic system explicitly in the first place. Much better would have been for your code to have simply written purely float coefficient data to data type=float[8] Matrix M and rhs b in the first place. Then pass those to LinearAlgebra:-LinearSolve. That can do a 2500x2500 linear system with a single rhs Vector in a few seconds if the data type is right.

@H-R You will not be able to solve the symbolic system for your FEM problem to get a general solution (which you would then try to instantiate at numeric values to get a solution for a particular problem). It's the wrong way to try it. A workable way is to produce a procedure which accepts numeric values for a particular problem, generates a numeric Matrix, and then solves that.

What do you intend on doing with the result? It will likely be enormous.

acer

The first assignment to parlist,

    [mu[p],seq(mu[p]+tau[p||j],j=3..K)]

does not depend on C and so could be pulled out of the vecC procedure. That procedure gets called C times, so it's a benefit to reduce what it does.

A few of the map calls in vecC can be rolled together. Each one creates temporary lists of expressions, and the less temporary expressions to garbage collect the better.

Some of the map calls could be done using only builtins and more than 2 parameters, rather than use of simple custom operators.

Right now vecC has a concatenation in eta[p||C] inside a mapped operator and gets done C times. How about making a local etapC of vecC which is assigned that right at the start of vecC and is used instead.

Lastly, are you taking a powerset of K things? The way you have it set up it first creates a list of multiplicands and then takes a powerset, and then combines those as products. Could you not take the call to `powerset` out of vecC, have it act on the initial list [mu[p],seq(mu[p]+tau[p||j],j=3..K)] of exponent addends, and then pass that to a reworked vecC which turns its input into the final products?

acer

@Jazen1 Are you supposed to use canned procedures from Maple which simply return the tangent and normal? Or are you supposed to show the steps of how it is done?

It may be that two square matrices can be so permuted only if they have the same Permanent. But is that sufficient as well as necessary? (proof might be in graph theoretic terms?)

And then there is generalization to nonsquare matrices to consider. Are your Matrices always square?

[edit: I think maybe this was a poor guess. I suspect that having the same permanent may be necessary but not sufficient.]

acer

@itsme I used a module in the proc I hit with Explore so that it could re-use the workspace Array for the many calls done while exploring. Even with the zeroing step on the Array it often wins to avoid producing Arrays for each call with must be memory managed as temp garbage. I could have made that re-usable workspace a global but a module keeps the global namespace cleaner.

The compiled procedure could be parallelized with Threads:-Task, following here, since the iterated map can be computed across entries of segments of the Array in a concurrent fashion. That's true even if each thread needs a few more iterations for the first entry in the segment. (A likely improvement to the iterated map proc would be to stop iterating, early, when the change in value for a given entry was smaller than a tolerance. I believe that my last revision already used the previous entrie's last value as the starting value for the next entrie's iteration.) The splitting by task would involve each thread computing on a subrange of the entire n=1..max_n.

So one could parallelize the inner workings of the iterated map code, acting inplace on an Array. To do this the Compiler:-Compile call would need to be passed the `threadsafe` option, otherwise the compiled call_external will be blocking. By blocking I mean that only one external call to the compiled external call is performed at any given time. The ability to Compile to non-blocking call_externals is a new feature of Maple 18, for which the external runtime of the Compiler is itself thread-safe.

So that's all about internally multithreading that inplace iterated mapping proc we'd been looking at. If you prefer to try and parallelize at a higher level then note that the resuable local workspace Array in an appliable module won't work (as I had it in the Explore thing). A local to a module is accessed globally by Threads. It is possible to create "thread-locals" inside a module but I consider that very advanced and have not tried it here. I'm talking here about parallelizing your unseen higher code which computes results for various values of the parameters (other than n). So without special coding that kind of parallelizing would need to create different result Arrays for each thread. The iterated mapping proc would still need to be Compiled with option `threadsafe`.

It's hard to say more usefully pertinent things, without knowing exactly what your full code is. It is often beneficial to highly optimize numeric stuff for serial computation, and quite often that can beat out parallelizing under-optimized code. That's not always true, but it pays to consider both approaches.