12

There are a couple of things that can be improved: Binary search can simply be replaced by a lookup table that is very cheap to build. Moreover, IntegerPart[i/len] leads to rational numbers (so no machine numbers), while Quotient can to the thing in vectorized way. Avoiding the use of a list of rules allows one to work with packed arrays. Using Pick instead ...


9

Using an actual SparseArray and machine precision this evaluates almost immediately with negligible memory footprint. s = 1; n = 10; b1 = N@Table[s - n, {n, 0, 2 s}]; σx = SparseArray@DiagonalMatrix[b1]; σI = IdentityMatrix[2 s + 1, SparseArray, WorkingPrecision -> MachinePrecision]; σxx[j_] := KroneckerProduct[## & @@ ConstantArray[σI, j - 1], σx, ## ...


7

AA=rand(3,3) tmp = spdiags(AA,0) tmp(2)=1 backmat=spdiags(tmp,0,AA); full(backmat) gives AA = 0.6948 0.0344 0.7655 0.3171 0.4387 0.7952 0.9502 0.3816 0.1869 tmp = 0.6948 0.4387 0.1869 tmp = 0.6948 1.0000 0.1869 ans = 0.6948 0.0344 0.7655 0.3171 1.0000 0.7952 0.9502 0.3816 ...


6

You don't need a library. This is already optimized in Mathematica. On my machine RepeatedTiming[t1 = H.H.H.H.H.H.H.H.H.H.f1;] RepeatedTiming[t2 = MatrixPower[H, 10, f1];] take 1.00s and 0.0032s respectively, for $>$ 300x speed-up. Note also that the slow part of your code is not the multiplication of a matrix and a vector, it's the multiplication of the ...


5

Not an answer, but some extended comments. Sparse matrix-sparse matrix multiplication MathematicaLover made the very good observation that--as written by OP--Mathematica performs sparse matrix-sparse matrix multiplies. In fact, Armadillo does it, too, it seems to be about twice as fast in that. I was quite surprised by the great margin by which Armadillo was ...


5

Here are two ways: dim = 7; lt = 2; s = SparseArray[{{i_, i_} :> 1.0 + I, Sequence @@ Table[{i, i + 1} -> 0.0, {i, lt}], Band[{2, 1}] -> 2.0, Band[{1, 2}] -> 2.0}, {dim, dim}, 0.0]; s // MatrixForm s2 = SparseArray[{{i_, i_} :> 1.0 + I, Band[{2, 1}] -> 2.0, Band[{lt + 1, lt + 2}] -> 2.0}, {dim, dim}, 0.0]; s == s2 (* ...


5

You seem to do something wrong because the LIL you provide is more suitable to assemble the transpose of the desired matrix in CRS format (or to assemble the desired matrix in CCS format). Since Mathematica uses CRS, I show you how to assemble the transpose. First two compiled helper functions: getColumnIndices = Compile[{{p, _Integer, 1}, {a, _Integer, 2}}, ...


4

Not an answer but too long for a comment. Well, this additional memory consumption is actually not as bad as you might think, and I am not so sure whether this should be considered a memory leak. The additive assembler has to sort the entries of Transpose[e[[1;;2]]]; to reorder the entries of e[[3]]; to compute the arrays for the row pointers; to compute ...


3

First OP's implementation with timing on my machine: First@AbsoluteTiming[ u1 = input; Do[f = SparseArray[{{i_} :> If[MemberQ[lll, i], 0, -2.0*u1[[i - 1]] - 2.0*u1[[i + 1]] - 2.0*u1[[i + dos]] - 2.0*u1[[i - dos]] + (1. + I)*u1[[i]]]}, {dos*dos}]; u1 = LinearSolve[s, f];, {j, 1, steps, 1}]; result0 = u1; ] 17.1247 The ...


2

Basically the following should do (replacing N by n and C by c): SparseArray[ Rule[ Transpose[{ Join @@ KeyValueMap[ConstantArray[#1, Length[#2]] &, n], Join @@ Values[n] }], Join @@ Values[c] ] ] To problem is that your numbers are too large. The matrix has as many rows as the largest key in the association n. And a sparse array has ...


2

Here are two versions of the same solution. First, using a general 2x2 matrix along the subdiagonal ClearAll[v] v[Nm_] := Block[{array, band}, array = With[{mat = Normal@SparseArray[Band[{2, 1}] -> Array[band, Nm - 1], Nm]}, ArrayFlatten[mat /. {0 -> ConstantArray[0, {2, 2}], band[n_] :> n{{b11, b12}, {b21, b22}}}]]; (array + ...


2

I don't know why the first case does not through an error, but SparseArray expectes the correct dimensions of the array as second argument. So V0 should read as follows: V0[Nm_] := SparseArray[{ Band[{3, 1}] -> Table[V[n], {n, 1, Nm - 1}], Band[{1, 3}] -> Table[ConjugateTranspose[V[n]], {n, 1, Nm - 1}] }, {2 Nm, 2 Nm} ]; The result ...


1

Error fixed in later version of Mathematica Given that others could not reproduce the error, it made sense that the problem was either in the version of Mathematica I was using or some very odd bug in my own computer. (Thank you to all who tried to reproduce the error.) Updating from 12.0 to 12.3 fixed the problem. It took a while to track down because it ...


1

The following (and my comment) overlooks the fact that u1 and f are dependent. Copy code, don't rewrite it -- also don't over edit once copied! When I look at Normal@f for various values of dos, I find the "edge" elements of f (considered as a dos by dos 2d grid) are 0 and the "inner" elements of f are 1.0I-7.0. Here, then, is a ...


Only top voted, non community-wiki answers of a minimum length are eligible