Skip to main content
added 3 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 6000;60000;
nc = 500;5000;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

118.480461406

0.00522008159

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 6000;
nc = 500;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

1.48046

0.00522

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 60000;
nc = 5000;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

18.1406

0.008159

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

added 1670 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 6000;
nc = 500;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

1.48046

0.00522

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 6000;
nc = 500;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

1.48046

0.00522

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

added 2 characters in body
Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.009902010061

Put as many operations as possible into one call:

mat[[ri]] += mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.009902

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Source Link
Henrik Schumacher
  • 109.4k
  • 7
  • 186
  • 322
Loading