|
|
@ -4,7 +4,7 @@ |
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
//
|
|
|
//
|
|
|
// Copyright (c) 2009-2010 Math.NET
|
|
|
// Copyright (c) 2009-2013 Math.NET
|
|
|
//
|
|
|
//
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
@ -28,9 +28,9 @@ |
|
|
// OTHER DEALINGS IN THE SOFTWARE.
|
|
|
// OTHER DEALINGS IN THE SOFTWARE.
|
|
|
// </copyright>
|
|
|
// </copyright>
|
|
|
|
|
|
|
|
|
|
|
|
using System; |
|
|
using MathNet.Numerics.Properties; |
|
|
using MathNet.Numerics.Properties; |
|
|
using MathNet.Numerics.Threading; |
|
|
using MathNet.Numerics.Threading; |
|
|
using System; |
|
|
|
|
|
|
|
|
|
|
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
{ |
|
|
{ |
|
|
@ -45,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|
|
/// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric
|
|
|
/// or positive definite, the constructor will throw an exception.
|
|
|
/// or positive definite, the constructor will throw an exception.
|
|
|
/// </remarks>
|
|
|
/// </remarks>
|
|
|
public class UserCholesky : Cholesky |
|
|
public sealed class UserCholesky : Cholesky |
|
|
{ |
|
|
{ |
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
/// Initializes a new instance of the <see cref="UserCholesky"/> class. This object will compute the
|
|
|
/// Initializes a new instance of the <see cref="UserCholesky"/> class. This object will compute the
|
|
|
@ -55,55 +55,57 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
|
|
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not a square matrix.</exception>
|
|
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
|
|
|
/// <exception cref="ArgumentException">If <paramref name="matrix"/> is not positive definite.</exception>
|
|
|
public UserCholesky(Matrix<Complex32> matrix) |
|
|
public static UserCholesky Create(Matrix<Complex32> matrix) |
|
|
{ |
|
|
{ |
|
|
if (matrix == null) |
|
|
|
|
|
{ |
|
|
|
|
|
throw new ArgumentNullException("matrix"); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
if (matrix.RowCount != matrix.ColumnCount) |
|
|
if (matrix.RowCount != matrix.ColumnCount) |
|
|
{ |
|
|
{ |
|
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|
|
throw new ArgumentException(Resources.ArgumentMatrixSquare); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
|
|
|
// Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
|
|
|
Factor = matrix.Clone(); |
|
|
var factor = matrix.Clone(); |
|
|
var tmpColumn = new Complex32[Factor.RowCount]; |
|
|
var tmpColumn = new Complex32[factor.RowCount]; |
|
|
|
|
|
|
|
|
// Main loop - along the diagonal
|
|
|
// Main loop - along the diagonal
|
|
|
for (var ij = 0; ij < Factor.RowCount; ij++) |
|
|
for (var ij = 0; ij < factor.RowCount; ij++) |
|
|
{ |
|
|
{ |
|
|
// "Pivot" element
|
|
|
// "Pivot" element
|
|
|
var tmpVal = Factor.At(ij, ij); |
|
|
var tmpVal = factor.At(ij, ij); |
|
|
|
|
|
|
|
|
if (tmpVal.Real > 0.0) |
|
|
if (tmpVal.Real > 0.0) |
|
|
{ |
|
|
{ |
|
|
tmpVal = tmpVal.SquareRoot(); |
|
|
tmpVal = tmpVal.SquareRoot(); |
|
|
Factor.At(ij, ij, tmpVal); |
|
|
factor.At(ij, ij, tmpVal); |
|
|
tmpColumn[ij] = tmpVal; |
|
|
tmpColumn[ij] = tmpVal; |
|
|
|
|
|
|
|
|
// Calculate multipliers and copy to local column
|
|
|
// Calculate multipliers and copy to local column
|
|
|
// Current column, below the diagonal
|
|
|
// Current column, below the diagonal
|
|
|
for (var i = ij + 1; i < Factor.RowCount; i++) |
|
|
for (var i = ij + 1; i < factor.RowCount; i++) |
|
|
{ |
|
|
{ |
|
|
Factor.At(i, ij, Factor.At(i, ij) / tmpVal); |
|
|
factor.At(i, ij, factor.At(i, ij)/tmpVal); |
|
|
tmpColumn[i] = Factor.At(i, ij); |
|
|
tmpColumn[i] = factor.At(i, ij); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
// Remaining columns, below the diagonal
|
|
|
// Remaining columns, below the diagonal
|
|
|
DoCholeskyStep(Factor, Factor.RowCount, ij + 1, Factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); |
|
|
DoCholeskyStep(factor, factor.RowCount, ij + 1, factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); |
|
|
} |
|
|
} |
|
|
else |
|
|
else |
|
|
{ |
|
|
{ |
|
|
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); |
|
|
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
for (var i = ij + 1; i < Factor.RowCount; i++) |
|
|
for (var i = ij + 1; i < factor.RowCount; i++) |
|
|
{ |
|
|
{ |
|
|
Factor.At(ij, i, Complex32.Zero); |
|
|
factor.At(ij, i, Complex32.Zero); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
return new UserCholesky(factor); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
UserCholesky(Matrix<Complex32> factor) |
|
|
|
|
|
: base(factor) |
|
|
|
|
|
{ |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
@ -115,14 +117,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
/// <param name="colLimit">Total columns</param>
|
|
|
/// <param name="colLimit">Total columns</param>
|
|
|
/// <param name="multipliers">Multipliers calculated previously</param>
|
|
|
/// <param name="multipliers">Multipliers calculated previously</param>
|
|
|
/// <param name="availableCores">Number of available processors</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) |
|
|
static void DoCholeskyStep(Matrix<Complex32> data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores) |
|
|
{ |
|
|
{ |
|
|
var tmpColCount = colLimit - firstCol; |
|
|
var tmpColCount = colLimit - firstCol; |
|
|
|
|
|
|
|
|
if ((availableCores > 1) && (tmpColCount > 200)) |
|
|
if ((availableCores > 1) && (tmpColCount > 200)) |
|
|
{ |
|
|
{ |
|
|
var tmpSplit = firstCol + (tmpColCount / 3); |
|
|
var tmpSplit = firstCol + (tmpColCount/3); |
|
|
var tmpCores = availableCores / 2; |
|
|
var tmpCores = availableCores/2; |
|
|
|
|
|
|
|
|
CommonParallel.Invoke( |
|
|
CommonParallel.Invoke( |
|
|
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), |
|
|
() => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), |
|
|
@ -135,7 +137,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
var tmpVal = multipliers[j]; |
|
|
var tmpVal = multipliers[j]; |
|
|
for (var i = j; i < rowDim; i++) |
|
|
for (var i = j; i < rowDim; i++) |
|
|
{ |
|
|
{ |
|
|
data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate())); |
|
|
data.At(i, j, data.At(i, j) - (multipliers[i]*tmpVal.Conjugate())); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
@ -148,17 +150,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|
|
public override void Solve(Matrix<Complex32> input, Matrix<Complex32> result) |
|
|
public override void Solve(Matrix<Complex32> input, Matrix<Complex32> result) |
|
|
{ |
|
|
{ |
|
|
if (input == null) |
|
|
|
|
|
{ |
|
|
|
|
|
throw new ArgumentNullException("input"); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
if (result == null) |
|
|
|
|
|
{ |
|
|
|
|
|
throw new ArgumentNullException("result"); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
// Check for proper dimensions.
|
|
|
|
|
|
if (result.RowCount != input.RowCount) |
|
|
if (result.RowCount != input.RowCount) |
|
|
{ |
|
|
{ |
|
|
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|
|
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); |
|
|
@ -186,10 +177,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
sum = result.At(i, c); |
|
|
sum = result.At(i, c); |
|
|
for (var k = i - 1; k >= 0; k--) |
|
|
for (var k = i - 1; k >= 0; k--) |
|
|
{ |
|
|
{ |
|
|
sum -= Factor.At(i, k) * result.At(k, c); |
|
|
sum -= Factor.At(i, k)*result.At(k, c); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
result.At(i, c, sum / Factor.At(i, i)); |
|
|
result.At(i, c, sum/Factor.At(i, i)); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
// Solve L'*X = Y;
|
|
|
// Solve L'*X = Y;
|
|
|
@ -198,10 +189,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
sum = result.At(i, c); |
|
|
sum = result.At(i, c); |
|
|
for (var k = i + 1; k < order; k++) |
|
|
for (var k = i + 1; k < order; k++) |
|
|
{ |
|
|
{ |
|
|
sum -= Factor.At(k, i).Conjugate() * result.At(k, c); |
|
|
sum -= Factor.At(k, i).Conjugate()*result.At(k, c); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
result.At(i, c, sum / Factor.At(i, i)); |
|
|
result.At(i, c, sum/Factor.At(i, i)); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
@ -213,18 +204,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|
|
public override void Solve(Vector<Complex32> input, Vector<Complex32> result) |
|
|
public override void Solve(Vector<Complex32> input, Vector<Complex32> result) |
|
|
{ |
|
|
{ |
|
|
// Check for proper arguments.
|
|
|
|
|
|
if (input == null) |
|
|
|
|
|
{ |
|
|
|
|
|
throw new ArgumentNullException("input"); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
if (result == null) |
|
|
|
|
|
{ |
|
|
|
|
|
throw new ArgumentNullException("result"); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
// Check for proper dimensions.
|
|
|
|
|
|
if (input.Count != result.Count) |
|
|
if (input.Count != result.Count) |
|
|
{ |
|
|
{ |
|
|
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|
|
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|
|
@ -245,10 +224,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
sum = result[i]; |
|
|
sum = result[i]; |
|
|
for (var k = i - 1; k >= 0; k--) |
|
|
for (var k = i - 1; k >= 0; k--) |
|
|
{ |
|
|
{ |
|
|
sum -= Factor.At(i, k) * result[k]; |
|
|
sum -= Factor.At(i, k)*result[k]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
result[i] = sum / Factor.At(i, i); |
|
|
result[i] = sum/Factor.At(i, i); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
// Solve L'*X = Y;
|
|
|
// Solve L'*X = Y;
|
|
|
@ -257,10 +236,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
sum = result[i]; |
|
|
sum = result[i]; |
|
|
for (var k = i + 1; k < order; k++) |
|
|
for (var k = i + 1; k < order; k++) |
|
|
{ |
|
|
{ |
|
|
sum -= Factor.At(k, i).Conjugate() * result[k]; |
|
|
sum -= Factor.At(k, i).Conjugate()*result[k]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
result[i] = sum / Factor.At(i, i); |
|
|
result[i] = sum/Factor.At(i, i); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
|
Christoph Ruegg
13 years ago