Skip to content

Commit f58711e

Browse files
cdrnetcdrnet
committed
Prefer looping over local arrays (II)
1 parent a56e095 commit f58711e

25 files changed

+563
-288
lines changed

‎src/Numerics/Integration/GaussLegendreRule.cs‎

Lines changed: 32 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -115,22 +115,25 @@ public static double Integrate(Func<double, double> f, double invervalBegin, dou
115115
double a = 0.5*(invervalEnd - invervalBegin);
116116
double b = 0.5*(invervalEnd + invervalBegin);
117117

118+
var weights = gaussLegendrePoint.Weights;
119+
var abscissas = gaussLegendrePoint.Abscissas;
120+
118121
if (order.IsOdd())
119122
{
120-
sum = gaussLegendrePoint.Weights[0]*f(b);
123+
sum = weights[0]*f(b);
121124
for (i = 1; i < m; i++)
122125
{
123-
ax = a*gaussLegendrePoint.Abscissas[i];
124-
sum += gaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
126+
ax = a*abscissas[i];
127+
sum += weights[i]*(f(b + ax) + f(b - ax));
125128
}
126129
}
127130
else
128131
{
129132
sum = 0.0;
130133
for (i = 0; i < m; i++)
131134
{
132-
ax = a*gaussLegendrePoint.Abscissas[i];
133-
sum += gaussLegendrePoint.Weights[i]*(f(b + ax) + f(b - ax));
135+
ax = a*abscissas[i];
136+
sum += weights[i]*(f(b + ax) + f(b - ax));
134137
}
135138
}
136139

@@ -157,22 +160,25 @@ public static Complex ContourIntegrate(Func<double, Complex> f, double invervalB
157160
double a = 0.5 * (invervalEnd - invervalBegin);
158161
double b = 0.5 * (invervalEnd + invervalBegin);
159162

163+
var weights = gaussLegendrePoint.Weights;
164+
var abscissas = gaussLegendrePoint.Abscissas;
165+
160166
if (order.IsOdd())
161167
{
162-
sum = gaussLegendrePoint.Weights[0] * f(b);
168+
sum = weights[0] * f(b);
163169
for (i = 1; i < m; i++)
164170
{
165-
ax = a * gaussLegendrePoint.Abscissas[i];
166-
sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax));
171+
ax = a * abscissas[i];
172+
sum += weights[i] * (f(b + ax) + f(b - ax));
167173
}
168174
}
169175
else
170176
{
171177
sum = 0.0;
172178
for (i = 0; i < m; i++)
173179
{
174-
ax = a * gaussLegendrePoint.Abscissas[i];
175-
sum += gaussLegendrePoint.Weights[i] * (f(b + ax) + f(b - ax));
180+
ax = a * abscissas[i];
181+
sum += weights[i] * (f(b + ax) + f(b - ax));
176182
}
177183
}
178184

@@ -202,34 +208,37 @@ public static double Integrate(Func<double, double, double> f, double invervalBe
202208
double c = 0.5*(invervalEndB - invervalBeginB);
203209
double d = 0.5*(invervalEndB + invervalBeginB);
204210

211+
var weights = gaussLegendrePoint.Weights;
212+
var abscissas = gaussLegendrePoint.Abscissas;
213+
205214
if (order.IsOdd())
206215
{
207-
sum = gaussLegendrePoint.Weights[0]*gaussLegendrePoint.Weights[0]*f(b, d);
216+
sum = weights[0]*weights[0]*f(b, d);
208217

209218
double t;
210219
for (j = 1, t = 0.0; j < m; j++)
211220
{
212-
cy = c*gaussLegendrePoint.Abscissas[j];
213-
t += gaussLegendrePoint.Weights[j]*(f(b, d + cy) + f(b, d - cy));
221+
cy = c*abscissas[j];
222+
t += weights[j]*(f(b, d + cy) + f(b, d - cy));
214223
}
215224

216-
sum += gaussLegendrePoint.Weights[0]*t;
225+
sum += weights[0]*t;
217226

218227
for (i = 1, t = 0.0; i < m; i++)
219228
{
220-
ax = a*gaussLegendrePoint.Abscissas[i];
221-
t += gaussLegendrePoint.Weights[i]*(f(b + ax, d) + f(b - ax, d));
229+
ax = a*abscissas[i];
230+
t += weights[i]*(f(b + ax, d) + f(b - ax, d));
222231
}
223232

224-
sum += gaussLegendrePoint.Weights[0]*t;
233+
sum += weights[0]*t;
225234

226235
for (i = 1; i < m; i++)
227236
{
228-
ax = a*gaussLegendrePoint.Abscissas[i];
237+
ax = a*abscissas[i];
229238
for (j = 1; j < m; j++)
230239
{
231-
cy = c*gaussLegendrePoint.Abscissas[j];
232-
sum += gaussLegendrePoint.Weights[i]*gaussLegendrePoint.Weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
240+
cy = c*abscissas[j];
241+
sum += weights[i]*weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
233242
}
234243
}
235244
}
@@ -238,11 +247,11 @@ public static double Integrate(Func<double, double, double> f, double invervalBe
238247
sum = 0.0;
239248
for (i = 0; i < m; i++)
240249
{
241-
ax = a*gaussLegendrePoint.Abscissas[i];
250+
ax = a*abscissas[i];
242251
for (j = 0; j < m; j++)
243252
{
244-
cy = c*gaussLegendrePoint.Abscissas[j];
245-
sum += gaussLegendrePoint.Weights[i]*gaussLegendrePoint.Weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
253+
cy = c*abscissas[j];
254+
sum += weights[i]*weights[j]*(f(b + ax, d + cy) + f(ax + b, d - cy) + f(b - ax, d + cy) + f(b - ax, d - cy));
246255
}
247256
}
248257
}

‎src/Numerics/Integration/GaussRule/GaussLegendrePointFactory.cs‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -89,17 +89,19 @@ static GaussPoint Map(double intervalBegin, double intervalEnd, GaussPoint gauss
8989

9090
int m = (gaussPoint.Order + 1) >> 1;
9191

92+
var gaussPointAbscissas = gaussPoint.Abscissas;
93+
var gaussPointWeights = gaussPoint.Weights;
9294
for (int i = 1; i <= m; i++)
9395
{
9496
int index1 = gaussPoint.Order - i;
9597
int index2 = i - 1;
9698
int index3 = m - i;
9799

98-
abscissas[index1] = gaussPoint.Abscissas[index3] * a + b;
99-
abscissas[index2] = -gaussPoint.Abscissas[index3] * a + b;
100+
abscissas[index1] = gaussPointAbscissas[index3] * a + b;
101+
abscissas[index2] = -gaussPointAbscissas[index3] * a + b;
100102

101-
weights[index1] = gaussPoint.Weights[index3] * a;
102-
weights[index2] = gaussPoint.Weights[index3] * a;
103+
weights[index1] = gaussPointWeights[index3] * a;
104+
weights[index2] = gaussPointWeights[index3] * a;
103105
}
104106

105107
return new GaussPoint(intervalBegin, intervalEnd, gaussPoint.Order, abscissas, weights);

‎src/Numerics/LinearAlgebra/Complex/DenseVector.cs‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -218,11 +218,12 @@ protected override void DoAdd(Complex scalar, Vector<Complex> result)
218218
{
219219
if (result is DenseVector dense)
220220
{
221+
var denseValues = dense._values;
221222
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
222223
{
223224
for (int i = a; i < b; i++)
224225
{
225-
dense._values[i] = _values[i] + scalar;
226+
denseValues[i] = _values[i] + scalar;
226227
}
227228
});
228229
}
@@ -276,11 +277,12 @@ protected override void DoSubtract(Complex scalar, Vector<Complex> result)
276277
{
277278
if (result is DenseVector dense)
278279
{
280+
var denseValues = dense._values;
279281
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
280282
{
281283
for (int i = a; i < b; i++)
282284
{
283-
dense._values[i] = _values[i] - scalar;
285+
denseValues[i] = _values[i] - scalar;
284286
}
285287
});
286288
}
@@ -412,10 +414,11 @@ protected override Complex DoConjugateDotProduct(Vector<Complex> other)
412414
{
413415
if (other is DenseVector denseVector)
414416
{
417+
var denseVectorValues = denseVector._values;
415418
var dot = Complex.Zero;
416419
for (var i = 0; i < _values.Length; i++)
417420
{
418-
dot += _values[i].Conjugate() * denseVector._values[i];
421+
dot += _values[i].Conjugate() * denseVectorValues[i];
419422
}
420423

421424
return dot;

‎src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -802,11 +802,12 @@ public override Matrix<Complex> Inverse()
802802
}
803803

804804
var inverse = (DiagonalMatrix)Clone();
805+
var inverseData = inverse._data;
805806
for (var i = 0; i < _data.Length; i++)
806807
{
807808
if (_data[i] != 0.0)
808809
{
809-
inverse._data[i] = 1.0 / _data[i];
810+
inverseData[i] = 1.0 / _data[i];
810811
}
811812
else
812813
{

‎src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs‎

Lines changed: 32 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -680,12 +680,15 @@ public override double InfinityNorm()
680680
/// <returns>The square root of the sum of the squared values.</returns>
681681
public override double FrobeniusNorm()
682682
{
683-
var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*ConjugateTranspose()).Storage;
684683
var norm = 0d;
684+
685+
var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*ConjugateTranspose()).Storage;
686+
var aatRowPointers = aat.RowPointers;
687+
var aatValues = aat.Values;
685688
for (var i = 0; i < aat.RowCount; i++)
686689
{
687-
var startIndex = aat.RowPointers[i];
688-
var endIndex = aat.RowPointers[i + 1];
690+
var startIndex = aatRowPointers[i];
691+
var endIndex = aatRowPointers[i + 1];
689692

690693
if (startIndex == endIndex)
691694
{
@@ -696,7 +699,7 @@ public override double FrobeniusNorm()
696699
{
697700
if (i == aat.ColumnIndices[j])
698701
{
699-
norm += aat.Values[j].Magnitude;
702+
norm += aatValues[j].Magnitude;
700703
}
701704
}
702705
}
@@ -742,13 +745,17 @@ protected override void DoAdd(Matrix<Complex> other, Matrix<Complex> result)
742745
}
743746

744747
var leftStorage = left._storage;
748+
var leftStorageRowPointers = leftStorage.RowPointers;
749+
var leftStorageColumnIndices = leftStorage.ColumnIndices;
750+
var leftStorageValues = leftStorage.Values;
751+
745752
for (var i = 0; i < leftStorage.RowCount; i++)
746753
{
747-
var endIndex = leftStorage.RowPointers[i + 1];
748-
for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
754+
var endIndex = leftStorageRowPointers[i + 1];
755+
for (var j = leftStorageRowPointers[i]; j < endIndex; j++)
749756
{
750-
var columnIndex = leftStorage.ColumnIndices[j];
751-
var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
757+
var columnIndex = leftStorageColumnIndices[j];
758+
var resVal = leftStorageValues[j] + result.At(i, columnIndex);
752759
result.At(i, columnIndex, resVal);
753760
}
754761
}
@@ -781,16 +788,19 @@ protected override void DoSubtract(Matrix<Complex> other, Matrix<Complex> result
781788
}
782789

783790
var otherStorage = sparseOther._storage;
791+
var otherStorageRowPointers = otherStorage.RowPointers;
792+
var otherStorageColumnIndices = otherStorage.ColumnIndices;
793+
var otherStorageValues = otherStorage.Values;
784794

785795
if (ReferenceEquals(this, sparseResult))
786796
{
787797
for (var i = 0; i < otherStorage.RowCount; i++)
788798
{
789-
var endIndex = otherStorage.RowPointers[i + 1];
790-
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
799+
var endIndex = otherStorageRowPointers[i + 1];
800+
for (var j = otherStorageRowPointers[i]; j < endIndex; j++)
791801
{
792-
var columnIndex = otherStorage.ColumnIndices[j];
793-
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
802+
var columnIndex = otherStorageColumnIndices[j];
803+
var resVal = sparseResult.At(i, columnIndex) - otherStorageValues[j];
794804
result.At(i, columnIndex, resVal);
795805
}
796806
}
@@ -913,6 +923,8 @@ protected override void DoMultiply(Matrix<Complex> other, Matrix<Complex> result
913923
var values = _storage.Values;
914924
if (other.Storage is DenseColumnMajorMatrixStorage<Complex> denseOther)
915925
{
926+
var denseOtherData = denseOther.Data;
927+
916928
// in this case we can directly address the underlying data-array
917929
for (var row = 0; row < RowCount; row++)
918930
{
@@ -930,7 +942,7 @@ protected override void DoMultiply(Matrix<Complex> other, Matrix<Complex> result
930942
var sum = Complex.Zero;
931943
for (var index = startIndex; index < endIndex; index++)
932944
{
933-
sum += values[index] * denseOther.Data[otherColumnStartPosition + columnIndices[index]];
945+
sum += values[index] * denseOtherData[otherColumnStartPosition + columnIndices[index]];
934946
}
935947

936948
result.At(row, column, sum);
@@ -1098,11 +1110,14 @@ protected override void DoTransposeAndMultiply(Matrix<Complex> other, Matrix<Com
10981110
var values = _storage.Values;
10991111

11001112
var otherStorage = otherSparse._storage;
1113+
var otherStorageRowPointers = otherStorage.RowPointers;
1114+
var otherStorageColumnIndices = otherStorage.ColumnIndices;
1115+
var otherStorageValues = otherStorage.Values;
11011116

11021117
for (var j = 0; j < RowCount; j++)
11031118
{
1104-
var startIndexOther = otherStorage.RowPointers[j];
1105-
var endIndexOther = otherStorage.RowPointers[j + 1];
1119+
var startIndexOther = otherStorageRowPointers[j];
1120+
var endIndexOther = otherStorageRowPointers[j + 1];
11061121

11071122
if (startIndexOther == endIndexOther)
11081123
{
@@ -1124,10 +1139,10 @@ protected override void DoTransposeAndMultiply(Matrix<Complex> other, Matrix<Com
11241139
var sum = Complex.Zero;
11251140
for (var index = startIndexOther; index < endIndexOther; index++)
11261141
{
1127-
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
1142+
var ind = _storage.FindItem(i, otherStorageColumnIndices[index]);
11281143
if (ind >= 0)
11291144
{
1130-
sum += otherStorage.Values[index] * values[ind];
1145+
sum += otherStorageValues[index] * values[ind];
11311146
}
11321147
}
11331148

0 commit comments

Comments
 (0)