# How to minmize the memory problem in a compile function for a given code?

n =3; (*Number of atoms*) L = 3 (*No=umber of atomic levels*);
CT[xk_]:=CT[xk]=ConjugateTranspose[xk]

Fx[xk_]:=Fx[xk]=FixedPoint[ArrayFlatten,xk]

gb[k_] = {{1},{0},{0}};
eb[k_] = {{0},{1},{0}};
rb[k_] = {{0}, {0}, {1}};
e1[l_] = Table[{KroneckerDelta[l, m]}, {m, L}];
σ[l_, m_] = e1[l].Transpose[e1[m]];

Hop[t_, k_] = g[t] gb[k].CT[eb[k]] + g[t] eb[k].CT[gb[k]]+ Ω[t] rb[k].CT[eb[k]] + Ω[t] eb[k].CT[rb[k]];

I1[i_Integer]:=I1[i]=IdentityMatrix[L^(i-1)];

I2[i_Integer]:=I2[i]=IdentityMatrix[L^(n-i)];

Hn[t_]=Sum[Fx[I1[i]\[TensorProduct]Hop[t,i]\[TensorProduct]I2[i]],{i,n}];

PrrρT[i_Integer]:= PrrρT[i]=Fx[I1[i]\[TensorProduct]Transpose[σ[2,3]]\[TensorProduct]I2[i]]

PeeρT[i_Integer]:= PeeρT[i]=Fx[I1[i]\[TensorProduct]Transpose[σ[1,2]]\[TensorProduct]I2[i]]

Prrρ[i_Integer]:= Prrρ[i]=Fx[I1[i]\[TensorProduct]σ[2,3]\[TensorProduct]I2[i]]

Peeρ[i_Integer]:= Peeρ[i]=Fx[I1[i]\[TensorProduct]σ[1,2]\[TensorProduct]I2[i]]

Hint[p_] = Sum[If[i < j, 10.0 (2 Ω0^2)/ Sqrt[2. Ω0^2 + Γe^2/4.]
PrrρT[i].Prrρ[i].PrrρT[j].Prrρ[j],0], {i, n}, {j, n}];

Hfinal[t_] = Hn[t] +  Hint[p];

Lfull[X_?MatrixQ,t_?NumericQ]=Sum[
(-Γr/2 X.#1.#2-Γr/2 #1.#2.X+Γr #2.X.#1
-Γe/2 X.#3.#4-Γe/2 #3.#4.X+Γe #4.X.#3
)&[PrrρT[i],Prrρ[i],PeeρT[i],Peeρ[i]]
,{i,n}];

(*Numerical Values*)
MHz = 2. N[π] 10^6; kHz = 2. N[π] 10^3;μs=1. 10^-6;

g[t_] = Ω0 N[Exp[-((t - tmax/2. - tmax/8.)^2/(2. (tmax/8)^2.))]];
Ω[t_] = g0 N[Exp[-((t - tmax/2. + tmax/8.)^2/(2. (tmax/8)^2.))]];
Γe = (38.0 MHz)/(2N[π]); Γr = (1. kHz)/(2N[π]);
Ω0 = 3. MHz; g0 = 3. MHz; tmax = 30. 10^-6;

(*Compile Function*)
m1 = With[{H1=Hfinal[t]},Compile[{{X, _Complex, 2}, {t, _Real}},
-I (H1.X - X.H1)+Lfull[X,t]]];

Fun[X_?MatrixQ,t_?NumericQ]:=m1[X,t];

mmu = MaxMemoryUsed[];

sol=NDSolve[{σ'[t]==Fun[σ[t],t],
σ[0]==Table[If[i == 1 ∧ j == 1,1.,0.],
{i,L^n},{j, L^n}]},σ[t],{t,0,tmax}, MaxStepSize -> 0.01,
MaxSteps -> 10^5 ];

mmu1=MaxMemoryUsed[]- mmu


If I consider n=6, then the NDSolve function is not working because of large memory problem...Please help me with this issue...

-
Why do you say "the compile function is not working"? For me it's the NDSolve that needs a lot of memory, so I stopped it to keep my computer operable. My impression is that you did try some optimizations in a quite arbitrary way, none of them seems to give any speedup, let alone safes you memory (which seems the main problem). Replacing all the Parallel functions with their sequential versions and using Compile with no options at all leaves the runtime exactly the same. You also might want to explain why you are using those options for NDSolve, can you explain why you chose them? – Albert Retey Jan 18 '14 at 19:18
@Albert, Thanks... I will try without using the parallel command. To clarify my statement... Compile function is working for n=3 case, but when I choose n=4 then it stops working and exit without compiling. You are right that NDSolve is consuming lots of memory. The method I am using for NDSolve is because otherwise it gives me an error and exit without solving for large n. & same is true for MaxSteps... – santosh Jan 18 '14 at 19:29
Which version/OS/compiler are you using? For me the Compile did work alright for n=4 with version 9.0.1 on Windows (both without any options and with your options). To make NDSolve finish without errors I only needed to increase MaxSteps, the rest of the options seems not to be necessary. Of course the quality of the solution might not be good enough without setting MaxStepSize, but using a small value for it will make the memory usage increase a lot (a factor of 3 for n=3). For n=4 I stopped NDSolve when the Kernel used more than 3 GB of memory... – Albert Retey Jan 18 '14 at 20:02
Thanks @Albert, without using parallel command its working for me as well and also Program is faster. But It is still not very efficient because I want to calculate for large n (like n=6). Somebody working on my field told me that they can get results in few hours for n=6 by using C programming. – santosh Jan 19 '14 at 19:50
My computer has 8 cores, but without using the parallel command, it only runs on one node. How I can Parallize my Program. – santosh Jan 20 '14 at 1:50

This is not really an answer but way too long for a comment. It should give some hints on how to tackle your problem and also show where the limits are...

# Speedup by Parallelizing:

To gain speedup from parallelism in most cases needs some fine-tuning: you only can get a 1/nprocs speedup in theory and in practice the communication overhead often eats up that speedup completely. That is especially true for a very high level language as Mathematica is, a naive parallelization will only give a speedup for very simple cases, which in a sense are "trivial" to parallelize. The code you did optimize is mainly operations on numeric matrices, which are already executed highly optimized by the underlying libraries and even will be automatically parallelized if appropriate, and in a much more efficient way than the high-level parallelism of the Parallel* functions. I wouldn't expect that you can gain much speedup from using the very high level parallelism here. It is more important to formulate the problem so that Mathematica and NDSolve can do their standard optimizations. I'm afraid that you won't really be able to make decent use of all the processors/cores on your machine for this problem, but you certainly will gain from every megabyte of RAM that it provides...

# Speedup by Compilation

If you compare your manually compiled version to one where m1 isn't compiled you will find that the compiled version is about twice as fast as the uncompiled, so compiling certainly helps. But NDSolve will usually automatically optimize and compile the equations before solving them if it can, and that will be even faster. For your code it can't compile because Fun is a black box for it due to the argument restrictions in its definition. If you instead use a very simple definition without any compilations you give NDSolve a chance to see what it really is that it's supposed to solve -- in this case that gives another factor two in speedup on my machine (and for much simple code :-)):

Fun[X_, t_] := -I (Hfinal[t].X - X.Hfinal[t]);


# NDSolve

For your specific code the problem is that cpu-time and memory are almost completely used in NDSolve and within it the possibilities to optimize are somewhat limited. On the other hand the code in NDSolve is already highly optimized and most critical parts of it already are compiled to machine code (most probably written in C). Nevertheless the generality of NDSolve has its price and it is very likely that the specialized C-code of your colleagues -- if written well -- will still outperform NDSolve by an order of magnitude or so. On top of that it often isn't really the computation time that gets problematic when solving large systems with NDSolve but the memory consumption and from what I tried with your code I'd expect the same here.

If I understand your problem correctly you are solving a relatively simple system of 3^n*3^n ODEs. For n=3 that are 729 equations but for n=6 that will be more than 500 thousend equations. That is certainly at the upper limit of what a system like Mathematica can reasonably be expected to solve, but it might be possible. The problem might not be the number of equations but that for some reason you need relatively small/many time steps and that makes it expensive, especially concerning memory.

One trick to make NDSolve use a lot less memory than by default is to not return the result as an interpolating function but only store those results you are actually interested in, e.g. not at every time step but only at certain times. The following call to NDSolve does that and stores the results only for 10 time steps:

mmu = MaxMemoryUsed[];
lastt = 0;
ClearAll@σsol;
Timing[
MemoryConstrained[
sol = NDSolve[{
σ'[t] == Fun[σ[t], t], σ[0] ==
Table[If[i == 1 && j == 1, 1., 0.], {i, L^n}, {j, L^n}]
}, {}, {t, 0, tmax}, MaxSteps -> 10^5,
StepMonitor :>
If[t - lastt > tmax/10, lastt = t; σsol[t] = σ[t];]
],
2^30(* ~1GB *)
];
]
mmu1 = MaxMemoryUsed[] - mmu


For me it runs only slightly slower than without the StepMonitor and for the n=3 case uses only 1MB additional memory compared to 500MB when returning the full interpolating function. With these optimizations I was able to run the n=4 case on my 2 years old laptop in about 5 minutes and with about 1GB memory used. I'd expect with that code you should be able to run the n=5 case on your machine.

The result is now stored as single matrices for the selected time steps in the downvalues of σsol. Here is an example how you can access and visualize the result:

tlist = DownValues[σsol][[All, 1, 1, 1]];
Manipulate[MatrixPlot[σsol[tlist[[k]]]], {k,1,Length[tlist],1}]


EDIT: as asked in a comment here is how you can safe only the elements on the diagonal, you just need to change the StepMonitor with the following:

 StepMonitor :> If[t - lastt > tmax/10, lastt = t; σsol[t] = Diagonal[σ[t]];]]


EDIT: To make a plot as described in your next comment you could then do this (are you really interested only in the sum of the first 5 entries of the diagonal?):

tlist = DownValues[σsol][[All, 1, 1, 1]];
ListLinePlot[{#,Total[Take[Re[σsol[#]], 5]]} & /@ tlist, PlotRange -> All]


To optimize even further I'd suggest to read tutorial/NDSolveOverview and see if you can fine tune the method-options. Using the methods and parameters that fit best instead of using the full magic that NDSolve uses to get a result could give some additional speedup. While it might be possible I wouldn't be surprised if it turns out that the n=6 case can't be solved with NDSolve even then on your machine. On the other hand, in some cases NDSolve is surprisingly performant. Whatever you try, it might turn out that using some specialized C-code to solve exactly that system can't be avoided and that that is written faster than finding the exact settings of NDSolve which would also work...

-
Thanks@Albert, I have an additional question. Is it possible to only evaluate diagonal elements of the matrix in my NDSolve? I only need to evaluate... p02:=Total[Table[Flatten[Re[σ[t]/.sol],1][[i,i]],{i,{1,2,4,5}}]]; p03:=p02+Total[Table[Flatten[Re[σ[t]/.sol],1][[i,i]],{i,{10,11,13,14}}]]; Plot[{p03,p13} ,{t,0,tmax},PlotRange->All] – santosh Jan 20 '14 at 17:53
@santosh: if I understand correctly that's easy enough, see my edit in the answer... – Albert Retey Jan 21 '14 at 21:33
Right, this is simple, but I don't want to plot Matrix. I want to plot this ... p02:=Total[Table[Flatten[Re[σ[t]/.sol],1][[i,i]],{i,{1,2,4,5}}]]; Plot[p03 ,{t,0,tmax},PlotRange->All] – santosh Jan 22 '14 at 0:09
In your case something like...Plot[σsol[tlist[[k]]], {k,1,Length[tlist]}] – santosh Jan 22 '14 at 0:56