Skip to content

Commit

Permalink
IRID-157: Matrix: Direct access to decomposition results (e.g. eigen …
Browse files Browse the repository at this point in the history
…values)

IRID-148: Vector: Scalar Addition/Subtraction Operators
Linear Algebra: double[] -> Vector
Minor Gendarme fixes.
  • Loading branch information
cdrnet committed May 28, 2008
1 parent e9bccd0 commit c4f8cf0
Show file tree
Hide file tree
Showing 7 changed files with 374 additions and 101 deletions.
Expand Up @@ -39,24 +39,24 @@ namespace MathNet.Numerics.LinearAlgebra
public class CholeskyDecomposition
{
/// <summary>Array for internal storage of decomposition.</summary>
double[][] L;
Matrix _l;

/// <summary>Symmetric and positive definite flag.</summary>
readonly bool _isSymmetricPositiveDefinite;

/// <summary>Cholesky algorithm for symmetric and positive definite matrix.</summary>
/// <param name="Arg">Square, symmetric matrix.</param>
/// <param name="m">Square, symmetric matrix.</param>
/// <returns>Structure to access L and isspd flag.</returns>
public
CholeskyDecomposition(
Matrix Arg
Matrix m
)
{
if(Arg.RowCount != Arg.ColumnCount)
if(m.RowCount != m.ColumnCount)
throw new InvalidOperationException(Resources.ArgumentMatrixSquare);

double[][] A = Arg;
L = Matrix.CreateMatrixData(Arg.RowCount, Arg.RowCount);
double[][] A = m;
double[][] L = Matrix.CreateMatrixData(m.RowCount, m.RowCount);

_isSymmetricPositiveDefinite = true; // ensure square

Expand Down Expand Up @@ -88,6 +88,8 @@ Matrix Arg
L[i][j] = 0.0;
}
}

_l = new Matrix(L);
}

/// <summary>Is the matrix symmetric and positive definite?</summary>
Expand All @@ -100,11 +102,20 @@ public bool IsSPD

/// <summary>Return triangular factor.</summary>
/// <returns>L</returns>
[Obsolete("Use the TriangularFactor property instead")]
public
Matrix
GetL()
{
return new Matrix(L);
return _l;
}

/// <summary>
/// Decomposition Triangular Factor Matrix (L).
/// </summary>
public Matrix TriangularFactor
{
get { return _l; }
}

/// <summary>Solve A*x = b</summary>
Expand All @@ -113,22 +124,24 @@ public bool IsSPD
/// <exception cref="System.ArgumentException">Matrix row dimensions must agree.</exception>
/// <exception cref="System.InvalidOperationException">Matrix is not symmetric positive definite.</exception>
public
double[]
Vector
Solve(
double[] b
Vector b
)
{
if(b.Length != L.Length)
if(b.Length != _l.RowCount)
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "b");
if(!_isSymmetricPositiveDefinite)
throw new InvalidOperationException(Resources.ArgumentMatrixSymmetricPositiveDefinite);

double[][] L = _l;
double[] bb = b;
double[] x = new double[b.Length];

// Solve L*y = b
for(int i = 0; i < x.Length; i++)
{
double sum = b[i];
double sum = bb[i];
for(int k = i - 1; k >= 0; k--)
{
sum -= L[i][k] * x[k];
Expand All @@ -147,7 +160,7 @@ public bool IsSPD
x[i] = sum / L[i][i];
}

return x;
return new Vector(x);
}

/// <summary>Solve A*X = B</summary>
Expand All @@ -161,12 +174,13 @@ public bool IsSPD
Matrix B
)
{
if(B.RowCount != L.Length)
if(B.RowCount != _l.RowCount)
throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension, "B");
if(!_isSymmetricPositiveDefinite)
throw new InvalidOperationException(Resources.ArgumentMatrixSymmetricPositiveDefinite);

int nx = B.ColumnCount;
double[][] L = _l;
double[][] BB = B;
double[][] X = Matrix.CreateMatrixData(L.Length, nx);

Expand Down
42 changes: 36 additions & 6 deletions src/app/MathNet.Iridium/Library/LinearAlgebra/ComplexVector.cs
Expand Up @@ -38,6 +38,7 @@ public class ComplexVector :
ICloneable
{

private Complex[] _data;
private int _length;

/// <summary>
Expand All @@ -48,11 +49,6 @@ public int Length
get { return _length; }
}

/// <summary>
/// Array for internal storage of elements.
/// </summary>
private Complex[] _data;

/// <summary>
/// Gets or sets the element indexed by <c>i</c>
/// in the <c>Vector</c>.
Expand Down Expand Up @@ -119,7 +115,7 @@ Complex value
/// the provided array as internal data structure.
/// </summary>
/// <param name="components">One-dimensional array of doubles.</param>
/// <seealso cref="Create"/>
/// <seealso cref="Create(Complex[])"/>
public
ComplexVector(
Complex[] components
Expand Down Expand Up @@ -150,6 +146,40 @@ Complex[] components
return new ComplexVector(newData);
}

/// <summary>
/// Constructs a complex vector from a real and an imaginary vector.
/// </summary>
/// <param name="realComponents">One-dimensional array of doubles representing the real part of the vector.</param>
/// <param name="imagComponents">One-dimensional array of doubles representing the imaginary part of the vector.</param>
public static
ComplexVector
Create(
IList<double> realComponents,
IList<double> imagComponents
)
{
if(null == realComponents)
{
throw new ArgumentNullException("realComponents");
}
if(null == imagComponents)
{
throw new ArgumentNullException("imagComponents");
}
if(realComponents.Count != imagComponents.Count)
{
throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths);
}

Complex[] newData = new Complex[realComponents.Count];
for(int i = 0; i < newData.Length; i++)
{
newData[i] = new Complex(realComponents[i], imagComponents[i]);
}

return new ComplexVector(newData);
}

// TODO: Add Random Vector Generations Methods here

/// <summary>
Expand Down

0 comments on commit c4f8cf0

Please sign in to comment.