# Efficient way for solving a matrix equation

I am looking for a more effcient method to perform the following analysis. I did my best, but running the following code on my machine takes 9.766s. I need to perform this analysis several times.

ClearAll["Global*"];

Jfunction =
Compile[{{NOunknownPar, _Integer, 0}, {dt, _Real,
0}, {Kgeneral, _Real, 2}, {Jacob1, _Real, 2}},
Module[{Ident, sj1, sj2, sj3, sf},
Ident = IdentityMatrix[2*NOunknownPar];
sj1 = Ident - (0.5*dt)*Kgeneral; sj2 = Ident + (0.5*dt)*Kgeneral;
sj3 = Dot @@ {sj2, Jacob1}; sf = LinearSolve[sj1, sj3]; sf],
CompilationTarget -> "C"];

NoVar = 79;
NoRep = 6001;
deltaT = 0.0001;
Jmat = IdentityMatrix[2*NoVar];
KK = ConstantArray[0, NoRep - 1];
Do[KK[[i]] = RandomReal[1, {2*NoVar, 2*NoVar}], {i, NoRep - 1}];

Do[Jmat = Jfunction[NoVar, deltaT, KK[[jJa]], Jmat], {jJa,
NoRep - 1}] // Timing

• There's not much use putting LinearSolve[] inside a compiled function. Also, sj3 = sj2.Jacob1 is much cleaner-looking. But, more importantly, how about stepping back from the code, and try to elaborate what you're trying to do here? – J. M. is away Jun 8 '15 at 9:53
• @Guesswhoitis Thanks for your comments. Probably it is not easy to explain exactly why I need this code, but generally I am trying to perform a non-linear dynamic analysis and in order to compute the Jacobian of the system with respect to the initial conditions, I need to solve a similar matrix equation at every time-step. – mak maak Jun 8 '15 at 10:05
• That last loop can be more cleanly done with FoldList[]. Also, KK = RandomReal[1, {NoRep - 1, 2*NoVar, 2*NoVar}]` is a better initialization. – J. M. is away Jun 8 '15 at 10:19