# Efficient ways to create matrices and solve matrix equations

I am attempting, for the first time, to use Mathematica to do some serious linear algebra. I would like to solve systems of equations of the form $$U_{n n'} f_{n'} = b_n.$$

I have an expression for $U_{n n'}$ that is of the general form $U_{n n'} = f_1(n) \delta_{n n'} + f_2(n,n')$. Is using 2 nested Table commands the easiest/most efficient way to build this matrix in Mathematica?

I am not entirely sure how large I will need to make the matrix (it results from discretizing an integral equation, so the number of rows/columns will be as many as I need to get an accurate solution). I guess that it could be as large as 10,000 x 10,000, maybe. Is LinearSolve efficient enough to handle these sized systems on a standard desktop PC? Is parallelization for this automatic or do I need to do something manually?

After I have found the solution, I am going to need to feed the solution to another equation to find the quantity that I am actually interested in. Is there anything I should do at the outset to make my life easier later?

I apologize for the general nature of my question, but this is all new ground for me, so I am not sure what general guidelines and practices are best.

• Depending on the number of elements that yield zero when operated on by f2 it might be useful to use SparseArray. Commented Apr 10, 2013 at 5:19
• Discretizing an integral equation will lead to very sparse matrices. You will want to therefore use SparseArray[] objects. Commented Apr 10, 2013 at 13:32
• @DanielLichtblau could you elaborate on why that would be the case? Based on the scheme I am using (approximating the integral as a finite sum, centered symmetrically about the origin) that matrices that I get are not at all sparse. Almost none of the entries are 0. Commented Apr 10, 2013 at 15:25
• Sorry, I was thinking "finite differencing". My mistake. Commented Apr 10, 2013 at 16:44

For making $U_{nn^\prime}$ you can also use the command Array (or ParallelArray if you want to utilize parallelization):

U = Array[f[1][#1] KroneckerDelta[#1,#2]+f[2][#1,#2]&,{n,m}]

LinearSolve is the safe bet, I dont know if it utilizes parallelization by default but you can always wrap it with Parallelize and check if Mathematica uses slave kernels. You can also give in LinearSolve just the matrix U and Mathematica will return a LinearSolveFunction object that you can apply to different vectors $b$ without the need to recalculate everything from the beginning.

What do you mean "Is there anything I should do at the outset to make my life easier later?" can you give an example?

The question is, why did you post a question, when you just could have tried it yourself.

n = 10^4;
m = RandomReal[{0, 1}, {n, n}];
b = RandomReal[{0, 1}, n];

LinearSolve[m, b]; // AbsoluteTiming

This needs 8 seconds here and the processor monitor suggests that it is at least partially parallelized

The limiting factor is surely your RAM. I have 32GB installed but n=5 stops already my kernel.

• Do you mean n=5 or n=10^5? Commented Apr 10, 2013 at 6:24
• The image of the process monitor was from the OS? If yes this is not the best way to look for parallelization because the process will probably change cpu or thread (context switching) and of course other processes would have been running on the same time. Parallel Kernel Status is a better way to monitor for parallelization. Commented Apr 10, 2013 at 6:56
• @Spawn1701D You are mistaken. The parallelisation which is used in built-in functions like matrix operations has absolutely nothing to do with Parallel Computing Tools of Mathematica. LinearSolve is executed on the main kernel and the parallelisation happens deep in the called mkl-library. Commented Apr 10, 2013 at 7:56
• @Spawn1701D Furthermore, the screenshot of my processor monitor was taken when nothing else was running. If I had taken a larger view you would have seen that in the time after LinearSolve finished there was almost no work load in all processors. Commented Apr 10, 2013 at 7:59
• @halirutan indeed but what you say is a low level multi threading performed utilizing intel cpus architecture I was referring to the high level parallelization feature of Mathematica. And, I assume, because it's just a multithreaded process you observe a partial parallelization. Commented Apr 10, 2013 at 8:33