ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

factor: merged Anrdiy's parallel Cholesky code and improved QR

la-knuth
Marcus Cuda 16 years ago
parent
0c969aff7c
commit
269f246509
9 changed files with 781 additions and 390 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 584
      src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
  2. 74
      src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
  3. 77
      src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
  4. 73
      src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs
  5. 77
      src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs
  6. 73
      src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs
  7. 71
      src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs
  8. 73
      src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs
  9. 69
      src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs

584
src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
View File

@ -229,7 +229,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
}
/// <summary>
/// Does a point wise division of two arrays <c>z = x / y</c>. This can be used
/// to divide elements of vectors or matrices.
@ -289,8 +288,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
case Norm.FrobeniusNorm:
break;
}
return ret;
throw new NotImplementedException();
}
/// <summary>
@ -1206,36 +1204,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("a");
}
for (var j = 0; j < order; j++)
var tmpColumn = new double[order];
// Main loop - along the diagonal
for (int ij = 0; ij < order; ij++)
{
var d = 0.0;
int index;
for (var k = 0; k < j; k++)
// "Pivot" element
double tmpVal = a[(ij * order) + ij];
if (tmpVal > 0.0)
{
var s = 0.0;
int i;
for (i = 0; i < k; i++)
tmpVal = Math.Sqrt(tmpVal);
a[(ij * order) + ij] = tmpVal;
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (int i = ij + 1; i < order; i++)
{
s += a[(i * order) + k] * a[(i * order) + j];
a[(ij * order) + i] /= tmpVal;
tmpColumn[i] = a[(ij * order) + i];
}
var tmp = k * order;
index = tmp + j;
a[index] = s = (a[index] - s) / a[tmp + k];
d += s * s;
// Remaining columns, below the diagonal
DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
index = (j * order) + j;
d = a[index] - d;
if (d <= 0.0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
a[index] = Math.Sqrt(d);
for (var k = j + 1; k < order; k++)
for (int i = ij + 1; i < order; i++)
{
a[(k * order) + j] = 0.0;
a[(i * order) + ij] = 0.0;
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliears calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(double[] data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data[(j * rowDim) + i] -= multipliers[i] * tmpVal;
}
}
}
}
@ -1428,13 +1464,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var minmn = Math.Min(rowsR, columnsR);
for (var i = 0; i < minmn; i++)
{
GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
GenerateColumn(work, r, rowsR, i, i);
ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
}
work[0] = rowsR * rowsR;
@ -1448,32 +1484,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="workIndex">Index of colunn in work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowCount">The number of rows</param>
/// <param name="rowStart">The first row in </param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowCount">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(double[] work, int workIndex, double[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnCount">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(double[] work, int workIndex, double[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
{
if (rowStart > rowEnd || columnStart > columnEnd)
if (rowStart > rowCount || columnStart > columnCount)
{
return;
}
var vector = new double[columnEnd - columnStart + 1];
for (var i = rowStart; i <= rowEnd; i++)
var tmpColCount = columnCount - columnStart;
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
() => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnCount; j++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * vector[j - columnStart];
var scale = 0.0;
for (var i = rowStart; i < rowCount; i++)
{
scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
for (var i = rowStart; i < rowCount; i++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * scale;
}
}
}
}
@ -1484,33 +1531,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="a">Initial matrix</param>
/// <param name="rowCount">The number of rows in matrix</param>
/// <param name="rowStart">The firts row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The firts row</param>
/// <param name="column">Column index</param>
private static void GenerateColumn(double[] work, double[] a, int rowCount, int rowStart, int rowEnd, int column)
private static void GenerateColumn(double[] work, double[] a, int rowCount, int row, int column)
{
var tmp = column * rowCount;
var index = tmp + rowStart;
var index = tmp + row;
CommonParallel.For(
rowStart,
rowEnd + 1,
row,
rowCount,
i =>
{
var iIndex = tmp + i;
work[iIndex - rowStart] = a[iIndex];
work[iIndex - row] = a[iIndex];
a[iIndex] = 0.0;
});
var norm = 0.0;
for (var i = 0; i < rowEnd - rowStart + 1; ++i)
for (var i = 0; i < rowCount - row; ++i)
{
var iindex = tmp + i;
norm += work[iindex] * work[iindex];
}
norm = Math.Sqrt(norm);
if (rowStart == rowEnd || norm == 0)
if (row == rowCount - 1 || norm == 0)
{
a[index] = -work[tmp];
work[tmp] = Math.Sqrt(2.0);
@ -1524,11 +1570,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
a[index] = -1.0 / scale;
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= scale);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale);
work[tmp] += 1.0;
var s = Math.Sqrt(1.0 / work[tmp]);
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= s);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s);
}
#endregion
@ -3977,43 +4023,81 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
public void CholeskyFactor(float[] a, int order)
{
var factor = new float[a.Length];
if (a == null)
{
throw new ArgumentNullException("a");
}
for (var j = 0; j < order; j++)
var tmpColumn = new float[order];
// Main loop - along the diagonal
for (var ij = 0; ij < order; ij++)
{
float d = 0.0f;
int index;
for (var k = 0; k < j; k++)
// "Pivot" element
var tmpVal = a[(ij * order) + ij];
if (tmpVal > 0.0)
{
float s = 0.0f;
int i;
for (i = 0; i < k; i++)
tmpVal = (float)Math.Sqrt(tmpVal);
a[(ij * order) + ij] = tmpVal;
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < order; i++)
{
s += factor[(i * order) + k] * factor[(i * order) + j];
a[(ij * order) + i] /= tmpVal;
tmpColumn[i] = a[(ij * order) + i];
}
var tmp = k * order;
index = tmp + j;
s = (a[index] - s) / factor[tmp + k];
factor[index] = s;
d += s * s;
// Remaining columns, below the diagonal
DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
index = (j * order) + j;
d = a[index] - d;
if (d <= 0.0F)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
factor[index] = (float)Math.Sqrt(d);
for (var k = j + 1; k < order; k++)
for (int i = ij + 1; i < order; i++)
{
factor[(k * order) + j] = 0.0F;
a[(i * order) + ij] = 0.0f;
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliears calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(float[] data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
Buffer.BlockCopy(factor, 0, a, 0, factor.Length * Constants.SizeOfFloat);
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data[(j * rowDim) + i] -= multipliers[i] * tmpVal;
}
}
}
}
/// <summary>
@ -4204,13 +4288,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var minmn = Math.Min(rowsR, columnsR);
for (var i = 0; i < minmn; i++)
{
GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
GenerateColumn(work, r, rowsR, i, i);
ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
}
work[0] = rowsR * rowsR;
@ -4224,32 +4308,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="workIndex">Index of colunn in work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowCount">The number of rows</param>
/// <param name="rowStart">The first row in </param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowCount">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(float[] work, int workIndex, float[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnCount">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(float[] work, int workIndex, float[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
{
if (rowStart > rowEnd || columnStart > columnEnd)
if (rowStart > rowCount || columnStart > columnCount)
{
return;
}
var vector = new float[columnEnd - columnStart + 1];
for (var i = rowStart; i <= rowEnd; i++)
var tmpColCount = columnCount - columnStart;
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
() => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnCount; j++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * vector[j - columnStart];
var scale = 0.0f;
for (var i = rowStart; i < rowCount; i++)
{
scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
for (var i = rowStart; i < rowCount; i++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * scale;
}
}
}
}
@ -4260,33 +4355,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="a">Initial matrix</param>
/// <param name="rowCount">The number of rows in matrix</param>
/// <param name="rowStart">The firts row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The firts row</param>
/// <param name="column">Column index</param>
private static void GenerateColumn(float[] work, float[] a, int rowCount, int rowStart, int rowEnd, int column)
private static void GenerateColumn(float[] work, float[] a, int rowCount, int row, int column)
{
var tmp = column * rowCount;
var index = tmp + rowStart;
var index = tmp + row;
CommonParallel.For(
rowStart,
rowEnd + 1,
row,
rowCount,
i =>
{
var iIndex = tmp + i;
work[iIndex - rowStart] = a[iIndex];
work[iIndex - row] = a[iIndex];
a[iIndex] = 0.0f;
});
var norm = 0.0;
for (var i = 0; i < rowEnd - rowStart + 1; ++i)
for (var i = 0; i < rowCount - row; ++i)
{
var iindex = tmp + i;
norm += work[iindex] * work[iindex];
}
norm = Math.Sqrt(norm);
if (rowStart == rowEnd || norm == 0)
if (row == rowCount - 1 || norm == 0)
{
a[index] = -work[tmp];
work[tmp] = (float)Math.Sqrt(2.0);
@ -4300,11 +4394,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
a[index] = -1.0f / scale;
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= scale);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale);
work[tmp] += 1.0f;
var s = (float)Math.Sqrt(1.0 / work[tmp]);
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= s);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s);
}
#endregion
@ -6769,35 +6863,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("a");
}
for (var i = 0; i < order; i++)
var tmpColumn = new Complex[order];
// Main loop - along the diagonal
for (var ij = 0; ij < order; ij++)
{
var d = Complex.Zero;
int index;
for (var j = 0; j < i; j++)
// "Pivot" element
var tmpVal = a[(ij * order) + ij];
if (tmpVal.Real > 0.0)
{
var s = Complex.Zero;
for (var k = 0; k < j; k++)
tmpVal = tmpVal.SquareRoot();
a[(ij * order) + ij] = tmpVal;
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < order; i++)
{
s += a[(k * order) + i] * a[(k * order) + j].Conjugate();
a[(ij * order) + i] /= tmpVal;
tmpColumn[i] = a[(ij * order) + i];
}
var tmp = j * order;
index = tmp + i;
a[index] = s = (a[index] - s) / a[tmp + j];
d += s * s.Conjugate();
// Remaining columns, below the diagonal
DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
index = (i * order) + i;
d = a[index] - d;
if (d.Real <= 0.0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
a[index] = d.SquareRoot();
for (var k = i + 1; k < order; k++)
for (var i = ij + 1; i < order; i++)
{
a[(k * order) + i] = 0.0;
a[(i * order) + ij] = 0.0;
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliears calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Complex[] data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data[(j * rowDim) + i] -= multipliers[i] * tmpVal.Conjugate();
}
}
}
}
@ -6990,13 +7123,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var minmn = Math.Min(rowsR, columnsR);
for (var i = 0; i < minmn; i++)
{
GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
GenerateColumn(work, r, rowsR, i, i);
ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
}
work[0] = rowsR * rowsR;
@ -7010,32 +7143,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="workIndex">Index of colunn in work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowCount">The number of rows</param>
/// <param name="rowStart">The first row in </param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowCount">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnCount">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
{
if (rowStart > rowEnd || columnStart > columnEnd)
if (rowStart > rowCount || columnStart > columnCount)
{
return;
}
var vector = new Complex[columnEnd - columnStart + 1];
for (var i = rowStart; i <= rowEnd; i++)
var tmpColCount = columnCount - columnStart;
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
() => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnCount; j++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * vector[j - columnStart];
var scale = Complex.Zero;
for (var i = rowStart; i < rowCount; i++)
{
scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
for (var i = rowStart; i < rowCount; i++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * scale;
}
}
}
}
@ -7046,33 +7190,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="a">Initial matrix</param>
/// <param name="rowCount">The number of rows in matrix</param>
/// <param name="rowStart">The firts row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The firts row</param>
/// <param name="column">Column index</param>
private static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int rowStart, int rowEnd, int column)
private static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int row, int column)
{
var tmp = column * rowCount;
var index = tmp + rowStart;
var index = tmp + row;
CommonParallel.For(
rowStart,
rowEnd + 1,
row,
rowCount,
i =>
{
var iIndex = tmp + i;
work[iIndex - rowStart] = a[iIndex];
work[iIndex - row] = a[iIndex];
a[iIndex] = Complex.Zero;
});
var norm = Complex.Zero;
for (var i = 0; i < rowEnd - rowStart + 1; ++i)
for (var i = 0; i < rowCount - row; ++i)
{
var index1 = tmp + i;
norm += work[index1].Magnitude * work[index1].Magnitude;
}
norm = norm.SquareRoot();
if (rowStart == rowEnd || norm.Magnitude == 0)
if (row == rowCount - 1 || norm.Magnitude == 0)
{
a[index] = -work[tmp];
work[tmp] = new Complex(2.0, 0).SquareRoot();
@ -7085,11 +7228,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
a[index] = -norm;
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] /= norm);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm);
work[tmp] += 1.0;
var s = (1.0 / work[tmp]).SquareRoot();
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
}
#endregion
@ -9505,35 +9648,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
throw new ArgumentNullException("a");
}
for (var i = 0; i < order; i++)
var tmpColumn = new Complex32[order];
// Main loop - along the diagonal
for (var ij = 0; ij < order; ij++)
{
var d = Complex32.Zero;
int index;
for (var j = 0; j < i; j++)
// "Pivot" element
var tmpVal = a[(ij * order) + ij];
if (tmpVal.Real > 0.0)
{
var s = Complex32.Zero;
for (var k = 0; k < j; k++)
tmpVal = tmpVal.SquareRoot();
a[(ij * order) + ij] = tmpVal;
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < order; i++)
{
s += a[(k * order) + i] * a[(k * order) + j].Conjugate();
a[(ij * order) + i] /= tmpVal;
tmpColumn[i] = a[(ij * order) + i];
}
var tmp = j * order;
index = tmp + i;
a[index] = s = (a[index] - s) / a[tmp + j];
d += s * s.Conjugate();
// Remaining columns, below the diagonal
DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
index = (i * order) + i;
d = a[index] - d;
if (d.Real <= 0.0f)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
a[index] = d.SquareRoot();
for (var k = i + 1; k < order; k++)
for (var i = ij + 1; i < order; i++)
{
a[(k * order) + i] = 0.0f;
a[(i * order) + ij] = 0.0f;
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliears calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Complex32[] data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data[(j * rowDim) + i] -= multipliers[i] * tmpVal.Conjugate();
}
}
}
}
@ -9726,13 +9908,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
var minmn = Math.Min(rowsR, columnsR);
for (var i = 0; i < minmn; i++)
{
GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
GenerateColumn(work, r, rowsR, i, i);
ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
}
work[0] = rowsR * rowsR;
@ -9746,32 +9928,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="workIndex">Index of colunn in work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowCount">The number of rows</param>
/// <param name="rowStart">The first row in </param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowCount">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(Complex32[] work, int workIndex, Complex32[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnCount">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(Complex32[] work, int workIndex, Complex32[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
{
if (rowStart > rowEnd || columnStart > columnEnd)
if (rowStart > rowCount || columnStart > columnCount)
{
return;
}
var vector = new Complex32[columnEnd - columnStart + 1];
for (var i = rowStart; i <= rowEnd; i++)
var tmpColCount = columnCount - columnStart;
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
() => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnCount; j++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * vector[j - columnStart];
var scale = Complex32.Zero;
for (var i = rowStart; i < rowCount; i++)
{
scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
}
for (var i = rowStart; i < rowCount; i++)
{
a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * scale;
}
}
}
}
@ -9782,33 +9975,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
/// <param name="work">Work array</param>
/// <param name="a">Initial matrix</param>
/// <param name="rowCount">The number of rows in matrix</param>
/// <param name="rowStart">The firts row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The firts row</param>
/// <param name="column">Column index</param>
private static void GenerateColumn(Complex32[] work, Complex32[] a, int rowCount, int rowStart, int rowEnd, int column)
private static void GenerateColumn(Complex32[] work, Complex32[] a, int rowCount, int row, int column)
{
var tmp = column * rowCount;
var index = tmp + rowStart;
var index = tmp + row;
CommonParallel.For(
rowStart,
rowEnd + 1,
row,
rowCount,
i =>
{
var iIndex = tmp + i;
work[iIndex - rowStart] = a[iIndex];
work[iIndex - row] = a[iIndex];
a[iIndex] = Complex32.Zero;
});
var norm = Complex32.Zero;
for (var i = 0; i < rowEnd - rowStart + 1; ++i)
for (var i = 0; i < rowCount - row; ++i)
{
var index1 = tmp + i;
norm += work[index1].Magnitude * work[index1].Magnitude;
}
norm = norm.SquareRoot();
if (rowStart == rowEnd || norm.Magnitude == 0)
if (row == rowCount - 1 || norm.Magnitude == 0)
{
a[index] = -work[tmp];
work[tmp] = new Complex32(2.0f, 0).SquareRoot();
@ -9821,11 +10013,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
}
a[index] = -norm;
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] /= norm);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm);
work[tmp] += 1.0f;
var s = (1.0f / work[tmp]).SquareRoot();
CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
}
#endregion

74
src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
View File

@ -33,8 +33,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
using System;
using System.Numerics;
using Generic;
using Generic.Factorization;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@ -69,32 +69,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
CholeskyFactor = matrix.Clone();
for (var i = 0; i < CholeskyFactor.RowCount; i++)
var tmpColumn = new Complex[CholeskyFactor.RowCount];
// Main loop - along the diagonal
for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
{
var d = Complex.Zero;
for (var j = 0; j < i; j++)
// "Pivot" element
var tmpVal = CholeskyFactor.At(ij, ij);
if (tmpVal.Real > 0.0)
{
var s = Complex.Zero;
for (var k = 0; k < j; k++)
tmpVal = tmpVal.SquareRoot();
CholeskyFactor.At(ij, ij, tmpVal);
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
s += CholeskyFactor.At(i, k) * CholeskyFactor.At(j, k).Conjugate();
CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
tmpColumn[i] = CholeskyFactor.At(i, ij);
}
s = (matrix.At(i, j) - s) / CholeskyFactor.At(j, j);
CholeskyFactor.At(i, j, s);
d += s * s.Conjugate();
// Remaining columns, below the diagonal
DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
d = matrix.At(i, i) - d;
if (d.Real <= 0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
CholeskyFactor.At(i, i, d.SquareRoot());
for (var k = i + 1; k < CholeskyFactor.RowCount; k++)
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
CholeskyFactor.At(i, k, 0.0);
CholeskyFactor.At(ij, i, Complex.Zero);
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliers calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Matrix<Complex> data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate()));
}
}
}
}

77
src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
View File

@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
using System.Numerics;
using Generic;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@ -77,13 +78,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
var u = new Complex[minmn][];
for (var i = 0; i < minmn; i++)
{
u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
u[i] = GenerateColumn(MatrixR, i, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
}
}
@ -91,27 +92,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
/// Generate column from initial matrix to work array
/// </summary>
/// <param name="a">Initial matrix</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The first row</param>
/// <param name="column">Column index</param>
/// <returns>Generated vector</returns>
private static Complex[] GenerateColumn(Matrix<Complex> a, int rowStart, int rowEnd, int column)
private static Complex[] GenerateColumn(Matrix<Complex> a, int row, int column)
{
var ru = rowEnd - rowStart + 1;
var ru = a.RowCount - row;
var u = new Complex[ru];
for (var i = rowStart; i <= rowEnd; i++)
for (var i = row; i < a.RowCount; i++)
{
u[i - rowStart] = a.At(i, column);
u[i - row] = a.At(i, column);
a.At(i, column, 0.0);
}
var norm = u.Aggregate(Complex.Zero, (current, t) => current + (t.Magnitude * t.Magnitude));
norm = norm.SquareRoot();
if (rowStart == rowEnd || norm.Magnitude == 0)
if (row == a.RowCount - 1 || norm.Magnitude == 0)
{
a.At(rowStart, column, -u[0]);
a.At(row, column, -u[0]);
u[0] = Math.Sqrt(2.0);
return u;
}
@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
norm = norm.Magnitude * (u[0] / u[0].Magnitude);
}
a.At(rowStart, column, -norm);
a.At(row, column, -norm);
for (var i = 0; i < ru; i++)
{
@ -140,40 +140,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
}
/// <summary>
/// Compute matrix Q or R
/// Perform calculation of Q or R
/// </summary>
/// <param name="u">H vectors values</param>
/// <param name="a">Matrix to calculate</param>
/// <param name="rowStart">Starting row index</param>
/// <param name="rowEnd">Ending row index</param>
/// <param name="columnStart">Starting column index</param>
/// <param name="columnEnd">Ending column index</param>
private static void ComputeQR(Complex[] u, Matrix<Complex> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="u">Work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowDim">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnDim">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(Complex[] u, Matrix<Complex> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
{
if (rowEnd < rowStart || columnEnd < columnStart)
if (rowDim < rowStart || columnDim < columnStart)
{
return;
}
var v = new Complex[columnEnd - columnStart + 1];
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = 0.0;
}
var tmpColCount = columnDim - columnStart;
for (var i = rowStart; i <= rowEnd; i++)
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
() => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnDim; j++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * v[j - columnStart]));
var scale = Complex.Zero;
for (var i = rowStart; i < rowDim; i++)
{
scale += u[i - rowStart] * a.At(i, j);
}
for (var i = rowStart; i < rowDim; i++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * scale));
}
}
}
}

73
src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs
View File

@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
using Generic;
using Numerics;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@ -68,32 +69,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
CholeskyFactor = matrix.Clone();
for (var i = 0; i < CholeskyFactor.RowCount; i++)
var tmpColumn = new Complex32[CholeskyFactor.RowCount];
// Main loop - along the diagonal
for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
{
var d = Complex32.Zero;
for (var j = 0; j < i; j++)
// "Pivot" element
var tmpVal = CholeskyFactor.At(ij, ij);
if (tmpVal.Real > 0.0)
{
var s = Complex32.Zero;
for (var k = 0; k < j; k++)
tmpVal = tmpVal.SquareRoot();
CholeskyFactor.At(ij, ij, tmpVal);
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
s += CholeskyFactor.At(i, k) * CholeskyFactor.At(j, k).Conjugate();
CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
tmpColumn[i] = CholeskyFactor.At(i, ij);
}
s = (matrix.At(i, j) - s) / CholeskyFactor.At(j, j);
CholeskyFactor.At(i, j, s);
d += s * s.Conjugate();
// Remaining columns, below the diagonal
DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
d = matrix.At(i, i) - d;
if (d.Real <= 0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
CholeskyFactor.At(i, i, d.SquareRoot());
for (var k = i + 1; k < CholeskyFactor.RowCount; k++)
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
CholeskyFactor.At(i, k, 0.0f);
CholeskyFactor.At(ij, i, Complex32.Zero);
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliers calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Matrix<Complex32> data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate()));
}
}
}
}

77
src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs
View File

@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
using Generic;
using Numerics;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@ -77,13 +78,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
var u = new Complex32[minmn][];
for (var i = 0; i < minmn; i++)
{
u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
u[i] = GenerateColumn(MatrixR, i, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
}
}
@ -91,27 +92,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
/// Generate column from initial matrix to work array
/// </summary>
/// <param name="a">Initial matrix</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The first row</param>
/// <param name="column">Column index</param>
/// <returns>Generated vector</returns>
private static Complex32[] GenerateColumn(Matrix<Complex32> a, int rowStart, int rowEnd, int column)
private static Complex32[] GenerateColumn(Matrix<Complex32> a, int row, int column)
{
var ru = rowEnd - rowStart + 1;
var ru = a.RowCount - row;
var u = new Complex32[ru];
for (var i = rowStart; i <= rowEnd; i++)
for (var i = row; i < a.RowCount; i++)
{
u[i - rowStart] = a.At(i, column);
u[i - row] = a.At(i, column);
a.At(i, column, 0.0f);
}
var norm = u.Aggregate(Complex32.Zero, (current, t) => current + (t.Magnitude * t.Magnitude));
norm = norm.SquareRoot();
if (rowStart == rowEnd || norm.Magnitude == 0)
if (row == a.RowCount - 1 || norm.Magnitude == 0)
{
a.At(rowStart, column, -u[0]);
a.At(row, column, -u[0]);
u[0] = (float)Math.Sqrt(2.0);
return u;
}
@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
norm = norm.Magnitude * (u[0] / u[0].Magnitude);
}
a.At(rowStart, column, -norm);
a.At(row, column, -norm);
for (var i = 0; i < ru; i++)
{
@ -140,40 +140,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
}
/// <summary>
/// Compute matrix Q or R
/// Perform calculation of Q or R
/// </summary>
/// <param name="u">H vectors values</param>
/// <param name="a">Matrix to calculate</param>
/// <param name="rowStart">Starting row index</param>
/// <param name="rowEnd">Ending row index</param>
/// <param name="columnStart">Starting column index</param>
/// <param name="columnEnd">Ending column index</param>
private static void ComputeQR(Complex32[] u, Matrix<Complex32> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="u">Work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowDim">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnDim">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(Complex32[] u, Matrix<Complex32> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
{
if (rowEnd < rowStart || columnEnd < columnStart)
if (rowDim < rowStart || columnDim < columnStart)
{
return;
}
var v = new Complex32[columnEnd - columnStart + 1];
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = 0.0f;
}
var tmpColCount = columnDim - columnStart;
for (var i = rowStart; i <= rowEnd; i++)
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
() => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnDim; j++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * v[j - columnStart]));
var scale = Complex32.Zero;
for (var i = rowStart; i < rowDim; i++)
{
scale += u[i - rowStart] * a.At(i, j);
}
for (var i = rowStart; i < rowDim; i++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * scale));
}
}
}
}

73
src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs
View File

@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
using System;
using Generic;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@ -67,32 +68,74 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
CholeskyFactor = matrix.Clone();
for (var j = 0; j < CholeskyFactor.RowCount; j++)
var tmpColumn = new double[CholeskyFactor.RowCount];
// Main loop - along the diagonal
for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
{
var d = 0.0;
for (var k = 0; k < j; k++)
// "Pivot" element
var tmpVal = CholeskyFactor.At(ij, ij);
if (tmpVal > 0.0)
{
var s = 0.0;
for (var i = 0; i < k; i++)
tmpVal = Math.Sqrt(tmpVal);
CholeskyFactor.At(ij, ij, tmpVal);
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
tmpColumn[i] = CholeskyFactor.At(i, ij);
}
s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
CholeskyFactor.At(j, k, s);
d += s * s;
// Remaining columns, below the diagonal
DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
d = matrix.At(j, j) - d;
if (d <= 0.0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
CholeskyFactor.At(j, j, Math.Sqrt(d));
for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
CholeskyFactor.At(j, k, 0.0);
CholeskyFactor.At(ij, i, 0.0);
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliers calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Matrix<double> data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal));
}
}
}
}

71
src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs
View File

@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
using System.Linq;
using Generic;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@ -76,41 +77,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
var u = new double[minmn][];
for (var i = 0; i < minmn; i++)
{
u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
u[i] = GenerateColumn(MatrixR, i, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
}
}
/// <summary>
/// Generate column from initial matrix to work array
/// </summary>
/// <param name="a">Initial matrix</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The firts row</param>
/// <param name="column">Column index</param>
/// <returns>Generated vector</returns>
private static double[] GenerateColumn(Matrix<double> a, int rowStart, int rowEnd, int column)
private static double[] GenerateColumn(Matrix<double> a, int row, int column)
{
var ru = rowEnd - rowStart + 1;
var ru = a.RowCount - row;
var u = new double[ru];
for (var i = rowStart; i <= rowEnd; i++)
for (var i = row; i < a.RowCount; i++)
{
u[i - rowStart] = a.At(i, rowStart);
a.At(i, rowStart, 0.0);
u[i - row] = a.At(i, row);
a.At(i, row, 0.0);
}
var norm = u.Sum(t => t * t);
norm = Math.Sqrt(norm);
if (rowStart == rowEnd || norm == 0)
if (row == a.RowCount - 1 || norm == 0)
{
a.At(rowStart, column, -u[0]);
a.At(row, column, -u[0]);
u[0] = Math.Sqrt(2.0);
return u;
}
@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
scale *= -1.0;
}
a.At(rowStart, column, -1.0 / scale);
a.At(row, column, -1.0 / scale);
for (var i = 0; i < ru; i++)
{
@ -145,35 +145,42 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
/// <param name="u">Work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowDim">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(double[] u, Matrix<double> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnDim">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(double[] u, Matrix<double> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
{
if (rowEnd < rowStart || columnEnd < columnStart)
if (rowDim < rowStart || columnDim < columnStart)
{
return;
}
var v = new double[columnEnd - columnStart + 1];
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = 0.0;
}
var tmpColCount = columnDim - columnStart;
for (var i = rowStart; i <= rowEnd; i++)
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
() => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnDim; j++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart]));
var scale = 0.0;
for (var i = rowStart; i < rowDim; i++)
{
scale += u[i - rowStart] * a.At(i, j);
}
for (var i = rowStart; i < rowDim; i++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart] * scale));
}
}
}
}

73
src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs
View File

@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
using System;
using Generic;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@ -67,32 +68,74 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
CholeskyFactor = matrix.Clone();
for (var j = 0; j < CholeskyFactor.RowCount; j++)
var tmpColumn = new float[CholeskyFactor.RowCount];
// Main loop - along the diagonal
for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
{
var d = 0.0;
for (var k = 0; k < j; k++)
// "Pivot" element
var tmpVal = CholeskyFactor.At(ij, ij);
if (tmpVal > 0.0)
{
var s = 0.0f;
for (var i = 0; i < k; i++)
tmpVal = (float)Math.Sqrt(tmpVal);
CholeskyFactor.At(ij, ij, tmpVal);
tmpColumn[ij] = tmpVal;
// Calculate multipliers and copy to local column
// Current column, below the diagonal
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
tmpColumn[i] = CholeskyFactor.At(i, ij);
}
s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
CholeskyFactor.At(j, k, s);
d += s * s;
// Remaining columns, below the diagonal
DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
}
d = matrix.At(j, j) - d;
if (d <= 0.0)
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
CholeskyFactor.At(j, j, (float)Math.Sqrt(d));
for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
{
CholeskyFactor.At(j, k, 0.0f);
CholeskyFactor.At(ij, i, 0.0f);
}
}
}
/// <summary>
/// Calculate Cholesky step
/// </summary>
/// <param name="data">Factor matrix</param>
/// <param name="rowDim">Number of rows</param>
/// <param name="firstCol">Column start</param>
/// <param name="colLimit">Total columns</param>
/// <param name="multipliers">Multipliers calculated previously</param>
/// <param name="availableCores">Number of available processors</param>
private static void DoCholeskyStep(Matrix<float> data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores)
{
var tmpColCount = colLimit - firstCol;
if ((availableCores > 1) && (tmpColCount > 200))
{
var tmpSplit = firstCol + (tmpColCount / 3);
var tmpCores = availableCores / 2;
CommonParallel.Invoke(
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
() => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
}
else
{
for (var j = firstCol; j < colLimit; j++)
{
var tmpVal = multipliers[j];
for (var i = j; i < rowDim; i++)
{
data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal));
}
}
}
}

69
src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs
View File

@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
using System.Linq;
using Generic;
using Properties;
using Threading;
/// <summary>
/// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@ -76,13 +77,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
var u = new float[minmn][];
for (var i = 0; i < minmn; i++)
{
u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
u[i] = GenerateColumn(MatrixR, i, i);
ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
}
}
@ -90,27 +91,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
/// Generate column from initial matrix to work array
/// </summary>
/// <param name="a">Initial matrix</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="row">The first row</param>
/// <param name="column">Column index</param>
/// <returns>Generated vector</returns>
private static float[] GenerateColumn(Matrix<float> a, int rowStart, int rowEnd, int column)
private static float[] GenerateColumn(Matrix<float> a, int row, int column)
{
var ru = rowEnd - rowStart + 1;
var ru = a.RowCount - row;
var u = new float[ru];
for (var i = rowStart; i <= rowEnd; i++)
for (var i = row; i < a.RowCount; i++)
{
u[i - rowStart] = a.At(i, rowStart);
a.At(i, rowStart, 0.0f);
u[i - row] = a.At(i, row);
a.At(i, row, 0.0f);
}
var norm = u.Sum(t => t * t);
norm = (float)Math.Sqrt(norm);
if (rowStart == rowEnd || norm == 0)
if (row == a.RowCount - 1 || norm == 0)
{
a.At(rowStart, column, -u[0]);
a.At(row, column, -u[0]);
u[0] = (float)Math.Sqrt(2.0);
return u;
}
@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
scale *= -1.0f;
}
a.At(rowStart, column, -1.0f / scale);
a.At(row, column, -1.0f / scale);
for (var i = 0; i < ru; i++)
{
@ -145,35 +145,42 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
/// <param name="u">Work array</param>
/// <param name="a">Q or R matrices</param>
/// <param name="rowStart">The first row</param>
/// <param name="rowEnd">The last row</param>
/// <param name="rowDim">The last row</param>
/// <param name="columnStart">The first column</param>
/// <param name="columnEnd">The last column</param>
private static void ComputeQR(float[] u, Matrix<float> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
/// <param name="columnDim">The last column</param>
/// <param name="availableCores">Number of available CPUs</param>
private static void ComputeQR(float[] u, Matrix<float> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
{
if (rowEnd < rowStart || columnEnd < columnStart)
if (rowDim < rowStart || columnDim < columnStart)
{
return;
}
var v = new float[columnEnd - columnStart + 1];
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = 0.0f;
}
var tmpColCount = columnDim - columnStart;
for (var i = rowStart; i <= rowEnd; i++)
if ((availableCores > 1) && (tmpColCount > 200))
{
for (var j = columnStart; j <= columnEnd; j++)
{
v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
}
}
var tmpSplit = columnStart + (tmpColCount / 2);
var tmpCores = availableCores / 2;
for (var i = rowStart; i <= rowEnd; i++)
CommonParallel.Invoke(
() => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
() => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
}
else
{
for (var j = columnStart; j <= columnEnd; j++)
for (var j = columnStart; j < columnDim; j++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart]));
var scale = 0.0f;
for (var i = rowStart; i < rowDim; i++)
{
scale += u[i - rowStart] * a.At(i, j);
}
for (var i = rowStart; i < rowDim; i++)
{
a.At(i, j, a.At(i, j) - (u[i - rowStart] * scale));
}
}
}
}

Write Preview
Loading…
Cancel
Save