// <copyright file="ManagedReferenceLinearAlgebraProvider.Complex32.cs" company="Math.NET">
|
// Math.NET Numerics, part of the Math.NET Project
|
// http://numerics.mathdotnet.com
|
// http://github.com/mathnet/mathnet-numerics
|
//
|
// Copyright (c) 2009-2020 Math.NET
|
//
|
// Permission is hereby granted, free of charge, to any person
|
// obtaining a copy of this software and associated documentation
|
// files (the "Software"), to deal in the Software without
|
// restriction, including without limitation the rights to use,
|
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
// copies of the Software, and to permit persons to whom the
|
// Software is furnished to do so, subject to the following
|
// conditions:
|
//
|
// The above copyright notice and this permission notice shall be
|
// included in all copies or substantial portions of the Software.
|
//
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
// OTHER DEALINGS IN THE SOFTWARE.
|
// </copyright>
|
|
using System;
|
using IStation.Numerics.Threading;
|
using Complex = System.Numerics.Complex;
|
using QRMethod = IStation.Numerics.LinearAlgebra.Factorization.QRMethod;
|
|
namespace IStation.Numerics.Providers.LinearAlgebra.ManagedReference
|
{
|
/// <summary>
|
/// The managed linear algebra provider.
|
/// </summary>
|
internal partial class ManagedReferenceLinearAlgebraProvider
|
{
|
/// <summary>
|
/// Adds a scaled vector to another: <c>result = y + alpha*x</c>.
|
/// </summary>
|
/// <param name="y">The vector to update.</param>
|
/// <param name="alpha">The value to scale <paramref name="x"/> by.</param>
|
/// <param name="x">The vector to add to <paramref name="y"/>.</param>
|
/// <param name="result">The result of the addition.</param>
|
/// <remarks>This is similar to the AXPY BLAS routine.</remarks>
|
public virtual void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (y.Length != x.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
if (alpha.IsZero())
|
{
|
y.Copy(result);
|
}
|
else if (alpha.IsOne())
|
{
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = y[i] + x[i];
|
}
|
});
|
}
|
else
|
{
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = y[i] + (alpha * x[i]);
|
}
|
});
|
}
|
|
}
|
|
/// <summary>
|
/// Scales an array. Can be used to scale a vector and a matrix.
|
/// </summary>
|
/// <param name="alpha">The scalar.</param>
|
/// <param name="x">The values to scale.</param>
|
/// <param name="result">This result of the scaling.</param>
|
/// <remarks>This is similar to the SCAL BLAS routine.</remarks>
|
public virtual void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result)
|
{
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (alpha.IsZero())
|
{
|
Array.Clear(result, 0, result.Length);
|
}
|
else if (alpha.IsOne())
|
{
|
x.Copy(result);
|
}
|
else
|
{
|
CommonParallel.For(0, x.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = alpha * x[i];
|
}
|
});
|
}
|
}
|
|
/// <summary>
|
/// Conjugates an array. Can be used to conjugate a vector and a matrix.
|
/// </summary>
|
/// <param name="x">The values to conjugate.</param>
|
/// <param name="result">This result of the conjugation.</param>
|
public virtual void ConjugateArray(Complex32[] x, Complex32[] result)
|
{
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
CommonParallel.For(0, x.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = x[i].Conjugate();
|
}
|
});
|
}
|
|
/// <summary>
|
/// Computes the dot product of x and y.
|
/// </summary>
|
/// <param name="x">The vector x.</param>
|
/// <param name="y">The vector y.</param>
|
/// <returns>The dot product of x and y.</returns>
|
/// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
|
public virtual Complex32 DotProduct(Complex32[] x, Complex32[] y)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (y.Length != x.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
var d = new Complex32(0.0F, 0.0F);
|
|
for (var i = 0; i < y.Length; i++)
|
{
|
d += y[i]*x[i];
|
}
|
|
return d;
|
}
|
|
/// <summary>
|
/// Does a point wise add of two arrays <c>z = x + y</c>. This can be used
|
/// to add vectors or matrices.
|
/// </summary>
|
/// <param name="x">The array x.</param>
|
/// <param name="y">The array y.</param>
|
/// <param name="result">The result of the addition.</param>
|
/// <remarks>There is no equivalent BLAS routine, but many libraries
|
/// provide optimized (parallel and/or vectorized) versions of this
|
/// routine.</remarks>
|
public virtual void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (y.Length != x.Length || y.Length != result.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = x[i] + y[i];
|
}
|
});
|
}
|
|
/// <summary>
|
/// Does a point wise subtraction of two arrays <c>z = x - y</c>. This can be used
|
/// to subtract vectors or matrices.
|
/// </summary>
|
/// <param name="x">The array x.</param>
|
/// <param name="y">The array y.</param>
|
/// <param name="result">The result of the subtraction.</param>
|
/// <remarks>There is no equivalent BLAS routine, but many libraries
|
/// provide optimized (parallel and/or vectorized) versions of this
|
/// routine.</remarks>
|
public virtual void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (y.Length != x.Length || y.Length != result.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = x[i] - y[i];
|
}
|
});
|
}
|
|
/// <summary>
|
/// Does a point wise multiplication of two arrays <c>z = x * y</c>. This can be used
|
/// to multiple elements of vectors or matrices.
|
/// </summary>
|
/// <param name="x">The array x.</param>
|
/// <param name="y">The array y.</param>
|
/// <param name="result">The result of the point wise multiplication.</param>
|
/// <remarks>There is no equivalent BLAS routine, but many libraries
|
/// provide optimized (parallel and/or vectorized) versions of this
|
/// routine.</remarks>
|
public virtual void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (y.Length != x.Length || y.Length != result.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = x[i] * y[i];
|
}
|
});
|
}
|
|
/// <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.
|
/// </summary>
|
/// <param name="x">The array x.</param>
|
/// <param name="y">The array y.</param>
|
/// <param name="result">The result of the point wise division.</param>
|
/// <remarks>There is no equivalent BLAS routine, but many libraries
|
/// provide optimized (parallel and/or vectorized) versions of this
|
/// routine.</remarks>
|
public virtual void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (y.Length != x.Length || y.Length != result.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = x[i] / y[i];
|
}
|
});
|
}
|
|
/// <summary>
|
/// Does a point wise power of two arrays <c>z = x ^ y</c>. This can be used
|
/// to raise elements of vectors or matrices to the powers of another vector or matrix.
|
/// </summary>
|
/// <param name="x">The array x.</param>
|
/// <param name="y">The array y.</param>
|
/// <param name="result">The result of the point wise power.</param>
|
/// <remarks>There is no equivalent BLAS routine, but many libraries
|
/// provide optimized (parallel and/or vectorized) versions of this
|
/// routine.</remarks>
|
public virtual void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result)
|
{
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (y.Length != x.Length || y.Length != result.Length)
|
{
|
throw new ArgumentException("All vectors must have the same dimensionality.");
|
}
|
|
CommonParallel.For(0, y.Length, 4096, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
result[i] = Complex32.Pow(x[i], y[i]);
|
}
|
});
|
}
|
|
/// <summary>
|
/// Computes the requested <see cref="Norm"/> of the matrix.
|
/// </summary>
|
/// <param name="norm">The type of norm to compute.</param>
|
/// <param name="rows">The number of rows.</param>
|
/// <param name="columns">The number of columns.</param>
|
/// <param name="matrix">The matrix to compute the norm from.</param>
|
/// <returns>The requested <see cref="Norm"/> of the matrix.</returns>
|
public virtual double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix)
|
{
|
switch (norm)
|
{
|
case Norm.OneNorm:
|
var norm1 = 0d;
|
for (var j = 0; j < columns; j++)
|
{
|
var s = 0d;
|
for (var i = 0; i < rows; i++)
|
{
|
s += matrix[(j*rows) + i].Magnitude;
|
}
|
|
norm1 = Math.Max(norm1, s);
|
}
|
return norm1;
|
case Norm.LargestAbsoluteValue:
|
var normMax = 0d;
|
for (var j = 0; j < columns; j++)
|
{
|
for (var i = 0; i < rows; i++)
|
{
|
normMax = Math.Max(matrix[(j * rows) + i].Magnitude, normMax);
|
}
|
}
|
return normMax;
|
case Norm.InfinityNorm:
|
var r = new double[rows];
|
for (var j = 0; j < columns; j++)
|
{
|
for (var i = 0; i < rows; i++)
|
{
|
r[i] += matrix[(j * rows) + i].Magnitude;
|
}
|
}
|
// TODO: reuse
|
var max = r[0];
|
for (int i = 0; i < r.Length; i++)
|
{
|
if (r[i] > max)
|
{
|
max = r[i];
|
}
|
}
|
return max;
|
case Norm.FrobeniusNorm:
|
var aat = new Complex32[rows*rows];
|
MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.ConjugateTranspose, 1.0f, matrix, rows, columns, matrix, rows, columns, 0.0f, aat);
|
var normF = 0d;
|
for (var i = 0; i < rows; i++)
|
{
|
normF += aat[(i * rows) + i].Magnitude;
|
}
|
return Math.Sqrt(normF);
|
default:
|
throw new NotSupportedException();
|
}
|
}
|
|
/// <summary>
|
/// Multiples two matrices. <c>result = x * y</c>
|
/// </summary>
|
/// <param name="x">The x matrix.</param>
|
/// <param name="rowsX">The number of rows in the x matrix.</param>
|
/// <param name="columnsX">The number of columns in the x matrix.</param>
|
/// <param name="y">The y matrix.</param>
|
/// <param name="rowsY">The number of rows in the y matrix.</param>
|
/// <param name="columnsY">The number of columns in the y matrix.</param>
|
/// <param name="result">Where to store the result of the multiplication.</param>
|
/// <remarks>This is a simplified version of the BLAS GEMM routine with alpha
|
/// set to 1.0 and beta set to 0.0, and x and y are not transposed.</remarks>
|
public virtual void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result)
|
{
|
// First check some basic requirement on the parameters of the matrix multiplication.
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (y == null)
|
{
|
throw new ArgumentNullException(nameof(y));
|
}
|
|
if (result == null)
|
{
|
throw new ArgumentNullException(nameof(result));
|
}
|
|
if (rowsX*columnsX != x.Length)
|
{
|
throw new ArgumentException("x.Length != xRows * xColumns");
|
}
|
|
if (rowsY*columnsY != y.Length)
|
{
|
throw new ArgumentException("y.Length != yRows * yColumns");
|
}
|
|
if (columnsX != rowsY)
|
{
|
throw new ArgumentException("xColumns != yRows");
|
}
|
|
if (rowsX*columnsY != result.Length)
|
{
|
throw new ArgumentException("xRows * yColumns != result.Length");
|
}
|
|
// Check whether we will be overwriting any of our inputs and make copies if necessary.
|
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
|
// as result, we can do it on a row wise basis. We should investigate this.
|
Complex32[] xdata;
|
if (ReferenceEquals(x, result))
|
{
|
xdata = (Complex32[]) x.Clone();
|
}
|
else
|
{
|
xdata = x;
|
}
|
|
Complex32[] ydata;
|
if (ReferenceEquals(y, result))
|
{
|
ydata = (Complex32[]) y.Clone();
|
}
|
else
|
{
|
ydata = y;
|
}
|
|
Array.Clear(result, 0, result.Length);
|
|
CacheObliviousMatrixMultiply(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, xdata, 0, 0, ydata, 0, 0, result, 0, 0, rowsX, columnsY, columnsX, rowsX, columnsY, columnsX, true);
|
}
|
|
/// <summary>
|
/// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
|
/// </summary>
|
/// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
|
/// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
|
/// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
|
/// <param name="a">The a matrix.</param>
|
/// <param name="rowsA">The number of rows in the <paramref name="a"/> matrix.</param>
|
/// <param name="columnsA">The number of columns in the <paramref name="a"/> matrix.</param>
|
/// <param name="b">The b matrix</param>
|
/// <param name="rowsB">The number of rows in the <paramref name="b"/> matrix.</param>
|
/// <param name="columnsB">The number of columns in the <paramref name="b"/> matrix.</param>
|
/// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
|
/// <param name="c">The c matrix.</param>
|
public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c)
|
{
|
int m; // The number of rows of matrix op(A) and of the matrix C.
|
int n; // The number of columns of matrix op(B) and of the matrix C.
|
int k; // The number of columns of matrix op(A) and the rows of the matrix op(B).
|
|
// First check some basic requirement on the parameters of the matrix multiplication.
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if ((int) transposeA > 111 && (int) transposeB > 111)
|
{
|
if (rowsA != columnsB)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
if (columnsA*rowsB != c.Length)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
m = columnsA;
|
n = rowsB;
|
k = rowsA;
|
}
|
else if ((int) transposeA > 111)
|
{
|
if (rowsA != rowsB)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
if (columnsA*columnsB != c.Length)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
m = columnsA;
|
n = columnsB;
|
k = rowsA;
|
}
|
else if ((int) transposeB > 111)
|
{
|
if (columnsA != columnsB)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
if (rowsA*rowsB != c.Length)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
m = rowsA;
|
n = rowsB;
|
k = columnsA;
|
}
|
else
|
{
|
if (columnsA != rowsB)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
if (rowsA*columnsB != c.Length)
|
{
|
throw new ArgumentOutOfRangeException();
|
}
|
|
m = rowsA;
|
n = columnsB;
|
k = columnsA;
|
}
|
|
if (alpha.IsZero() && beta.IsZero())
|
{
|
Array.Clear(c, 0, c.Length);
|
return;
|
}
|
|
// Check whether we will be overwriting any of our inputs and make copies if necessary.
|
// TODO - we can don't have to allocate a completely new matrix when x or y point to the same memory
|
// as result, we can do it on a row wise basis. We should investigate this.
|
Complex32[] adata;
|
if (ReferenceEquals(a, c))
|
{
|
adata = (Complex32[]) a.Clone();
|
}
|
else
|
{
|
adata = a;
|
}
|
|
Complex32[] bdata;
|
if (ReferenceEquals(b, c))
|
{
|
bdata = (Complex32[]) b.Clone();
|
}
|
else
|
{
|
bdata = b;
|
}
|
|
if (beta.IsZero())
|
{
|
Array.Clear(c, 0, c.Length);
|
}
|
else if (!beta.IsOne())
|
{
|
ScaleArray(beta, c, c);
|
}
|
|
if (alpha.IsZero())
|
{
|
return;
|
}
|
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, adata, 0, 0, bdata, 0, 0, c, 0, 0, m, n, k, m, n, k, true);
|
}
|
|
/// <summary>
|
/// Cache-Oblivious Matrix Multiplication
|
/// </summary>
|
/// <param name="transposeA">if set to <c>true</c> transpose matrix A.</param>
|
/// <param name="transposeB">if set to <c>true</c> transpose matrix B.</param>
|
/// <param name="alpha">The value to scale the matrix A with.</param>
|
/// <param name="matrixA">The matrix A.</param>
|
/// <param name="shiftArow">Row-shift of the left matrix</param>
|
/// <param name="shiftAcol">Column-shift of the left matrix</param>
|
/// <param name="matrixB">The matrix B.</param>
|
/// <param name="shiftBrow">Row-shift of the right matrix</param>
|
/// <param name="shiftBcol">Column-shift of the right matrix</param>
|
/// <param name="result">The matrix C.</param>
|
/// <param name="shiftCrow">Row-shift of the result matrix</param>
|
/// <param name="shiftCcol">Column-shift of the result matrix</param>
|
/// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
|
/// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
|
/// <param name="k">The number of columns of matrix op(A) and the rows of the matrix op(B).</param>
|
/// <param name="constM">The constant number of rows of matrix op(A) and of the matrix C.</param>
|
/// <param name="constN">The constant number of columns of matrix op(B) and of the matrix C.</param>
|
/// <param name="constK">The constant number of columns of matrix op(A) and the rows of the matrix op(B).</param>
|
/// <param name="first">Indicates if this is the first recursion.</param>
|
static void CacheObliviousMatrixMultiply(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] matrixA, int shiftArow, int shiftAcol, Complex32[] matrixB, int shiftBrow, int shiftBcol, Complex32[] result, int shiftCrow, int shiftCcol, int m, int n, int k, int constM, int constN, int constK, bool first)
|
{
|
if (m + n <= Control.ParallelizeOrder || m == 1 || n == 1 || k == 1)
|
{
|
if ((int) transposeA > 111 && (int) transposeB > 111)
|
{
|
if ((int) transposeA > 112 && (int) transposeB > 112)
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
else if ((int) transposeA > 112)
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
else if ((int) transposeB > 112)
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
else
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
}
|
else if ((int) transposeA > 111)
|
{
|
if ((int) transposeA > 112)
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol].Conjugate()*
|
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
else
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[(matArowPos*constK) + k1 + shiftAcol]*
|
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
}
|
else if ((int) transposeB > 111)
|
{
|
if ((int) transposeB > 112)
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos].Conjugate();
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
else
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
|
matrixB[((k1 + shiftBrow)*constN) + matBcolPos];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
}
|
else
|
{
|
for (var m1 = 0; m1 < m; m1++)
|
{
|
var matArowPos = m1 + shiftArow;
|
var matCrowPos = m1 + shiftCrow;
|
for (var n1 = 0; n1 < n; ++n1)
|
{
|
var matBcolPos = n1 + shiftBcol;
|
var sum = Complex32.Zero;
|
for (var k1 = 0; k1 < k; ++k1)
|
{
|
sum += matrixA[((k1 + shiftAcol)*constM) + matArowPos]*
|
matrixB[(matBcolPos*constK) + k1 + shiftBrow];
|
}
|
|
result[((n1 + shiftCcol)*constM) + matCrowPos] += alpha*sum;
|
}
|
}
|
}
|
}
|
else
|
{
|
// divide and conquer
|
int m2 = m/2, n2 = n/2, k2 = k/2;
|
|
if (first)
|
{
|
CommonParallel.Invoke(
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false));
|
|
CommonParallel.Invoke(
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false),
|
() => CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false));
|
}
|
else
|
{
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k2, constM, constN, constK, false);
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k2, constM, constN, constK, false);
|
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow, shiftCcol, m2, n2, k - k2, constM, constN, constK, false);
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow, shiftCcol + n2, m2, n - n2, k - k2, constM, constN, constK, false);
|
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k2, constM, constN, constK, false);
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol, matrixB, shiftBrow, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k2, constM, constN, constK, false);
|
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol, result, shiftCrow + m2, shiftCcol, m - m2, n2, k - k2, constM, constN, constK, false);
|
CacheObliviousMatrixMultiply(transposeA, transposeB, alpha, matrixA, shiftArow + m2, shiftAcol + k2, matrixB, shiftBrow + k2, shiftBcol + n2, result, shiftCrow + m2, shiftCcol + n2, m - m2, n - n2, k - k2, constM, constN, constK, false);
|
}
|
}
|
}
|
|
/// <summary>
|
/// Computes the LUP factorization of A. P*A = L*U.
|
/// </summary>
|
/// <param name="data">An <paramref name="order"/> by <paramref name="order"/> matrix. The matrix is overwritten with the
|
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of <paramref name="data"/> (the diagonal is always 1.0
|
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of <paramref name="data"/>.</param>
|
/// <param name="order">The order of the square matrix <paramref name="data"/>.</param>
|
/// <param name="ipiv">On exit, it contains the pivot indices. The size of the array must be <paramref name="order"/>.</param>
|
/// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
|
public virtual void LUFactor(Complex32[] data, int order, int[] ipiv)
|
{
|
if (data == null)
|
{
|
throw new ArgumentNullException(nameof(data));
|
}
|
|
if (ipiv == null)
|
{
|
throw new ArgumentNullException(nameof(ipiv));
|
}
|
|
if (data.Length != order*order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(data));
|
}
|
|
if (ipiv.Length != order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(ipiv));
|
}
|
|
// Initialize the pivot matrix to the identity permutation.
|
for (var i = 0; i < order; i++)
|
{
|
ipiv[i] = i;
|
}
|
|
var vecLUcolj = new Complex32[order];
|
|
// Outer loop.
|
for (var j = 0; j < order; j++)
|
{
|
var indexj = j*order;
|
var indexjj = indexj + j;
|
|
// Make a copy of the j-th column to localize references.
|
for (var i = 0; i < order; i++)
|
{
|
vecLUcolj[i] = data[indexj + i];
|
}
|
|
// Apply previous transformations.
|
for (var i = 0; i < order; i++)
|
{
|
// Most of the time is spent in the following dot product.
|
var kmax = Math.Min(i, j);
|
var s = Complex32.Zero;
|
for (var k = 0; k < kmax; k++)
|
{
|
s += data[(k*order) + i]*vecLUcolj[k];
|
}
|
|
data[indexj + i] = vecLUcolj[i] -= s;
|
}
|
|
// Find pivot and exchange if necessary.
|
var p = j;
|
for (var i = j + 1; i < order; i++)
|
{
|
if (vecLUcolj[i].Magnitude > vecLUcolj[p].Magnitude)
|
{
|
p = i;
|
}
|
}
|
|
if (p != j)
|
{
|
for (var k = 0; k < order; k++)
|
{
|
var indexk = k*order;
|
var indexkp = indexk + p;
|
var indexkj = indexk + j;
|
var temp = data[indexkp];
|
data[indexkp] = data[indexkj];
|
data[indexkj] = temp;
|
}
|
|
ipiv[j] = p;
|
}
|
|
// Compute multipliers.
|
if (j < order & data[indexjj] != 0.0f)
|
{
|
for (var i = j + 1; i < order; i++)
|
{
|
data[indexj + i] /= data[indexjj];
|
}
|
}
|
}
|
}
|
|
/// <summary>
|
/// Computes the inverse of matrix using LU factorization.
|
/// </summary>
|
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
|
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
|
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
|
public virtual void LUInverse(Complex32[] a, int order)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (a.Length != order*order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(a));
|
}
|
|
var ipiv = new int[order];
|
LUFactor(a, order, ipiv);
|
LUInverseFactored(a, order, ipiv);
|
}
|
|
/// <summary>
|
/// Computes the inverse of a previously factored matrix.
|
/// </summary>
|
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
|
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
|
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
|
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
|
public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (ipiv == null)
|
{
|
throw new ArgumentNullException(nameof(ipiv));
|
}
|
|
if (a.Length != order*order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(a));
|
}
|
|
if (ipiv.Length != order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(ipiv));
|
}
|
|
var inverse = new Complex32[a.Length];
|
for (var i = 0; i < order; i++)
|
{
|
inverse[i + (order*i)] = Complex32.One;
|
}
|
|
LUSolveFactored(order, a, order, ipiv, inverse);
|
inverse.Copy(a);
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using LU factorization.
|
/// </summary>
|
/// <param name="columnsOfB">The number of columns of B.</param>
|
/// <param name="a">The square matrix A.</param>
|
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
|
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
|
/// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
|
public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (a.Length != order*order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(a));
|
}
|
|
if (b.Length != order*columnsOfB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (ReferenceEquals(a, b))
|
{
|
throw new ArgumentException("Arguments must be different objects.");
|
}
|
|
var ipiv = new int[order];
|
var clone = new Complex32[a.Length];
|
a.Copy(clone);
|
LUFactor(clone, order, ipiv);
|
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using a previously factored A matrix.
|
/// </summary>
|
/// <param name="columnsOfB">The number of columns of B.</param>
|
/// <param name="a">The factored A matrix.</param>
|
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
|
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
|
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
|
/// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
|
public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (ipiv == null)
|
{
|
throw new ArgumentNullException(nameof(ipiv));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (a.Length != order*order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(a));
|
}
|
|
if (ipiv.Length != order)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(ipiv));
|
}
|
|
if (b.Length != order*columnsOfB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (ReferenceEquals(a, b))
|
{
|
throw new ArgumentException("Arguments must be different objects.");
|
}
|
|
// Compute the column vector P*B
|
for (var i = 0; i < ipiv.Length; i++)
|
{
|
if (ipiv[i] == i)
|
{
|
continue;
|
}
|
|
var p = ipiv[i];
|
for (var j = 0; j < columnsOfB; j++)
|
{
|
var indexk = j*order;
|
var indexkp = indexk + p;
|
var indexkj = indexk + i;
|
var temp = b[indexkp];
|
b[indexkp] = b[indexkj];
|
b[indexkj] = temp;
|
}
|
}
|
|
// Solve L*Y = P*B
|
for (var k = 0; k < order; k++)
|
{
|
var korder = k*order;
|
for (var i = k + 1; i < order; i++)
|
{
|
for (var j = 0; j < columnsOfB; j++)
|
{
|
var index = j*order;
|
b[i + index] -= b[k + index]*a[i + korder];
|
}
|
}
|
}
|
|
// Solve U*X = Y;
|
for (var k = order - 1; k >= 0; k--)
|
{
|
var korder = k + (k*order);
|
for (var j = 0; j < columnsOfB; j++)
|
{
|
b[k + (j*order)] /= a[korder];
|
}
|
|
korder = k*order;
|
for (var i = 0; i < k; i++)
|
{
|
for (var j = 0; j < columnsOfB; j++)
|
{
|
var index = j*order;
|
b[i + index] -= b[k + index]*a[i + korder];
|
}
|
}
|
}
|
}
|
|
/// <summary>
|
/// Computes the Cholesky factorization of A.
|
/// </summary>
|
/// <param name="a">On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
|
/// the Cholesky factorization.</param>
|
/// <param name="order">The number of rows or columns in the matrix.</param>
|
/// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
|
public virtual void CholeskyFactor(Complex32[] a, int order)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
var tmpColumn = new Complex32[order];
|
|
// Main loop - along the diagonal
|
for (var ij = 0; ij < order; ij++)
|
{
|
// "Pivot" element
|
var tmpVal = a[(ij*order) + ij];
|
|
if (tmpVal.Real > 0.0)
|
{
|
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++)
|
{
|
a[(ij*order) + i] /= tmpVal;
|
tmpColumn[i] = a[(ij*order) + i];
|
}
|
|
// Remaining columns, below the diagonal
|
DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.MaxDegreeOfParallelism);
|
}
|
else
|
{
|
throw new ArgumentException("Matrix must be positive definite.");
|
}
|
|
for (var i = ij + 1; i < order; i++)
|
{
|
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">Multipliers calculated previously</param>
|
/// <param name="availableCores">Number of available processors</param>
|
static void DoCholeskyStep(Complex32[] data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
|
{
|
var tmpColCount = colLimit - firstCol;
|
|
if ((availableCores > 1) && (tmpColCount > Control.ParallelizeElements))
|
{
|
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();
|
}
|
}
|
}
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using Cholesky factorization.
|
/// </summary>
|
/// <param name="a">The square, positive definite matrix A.</param>
|
/// <param name="orderA">The number of rows and columns in A.</param>
|
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
|
/// <param name="columnsB">The number of columns in the B matrix.</param>
|
/// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
|
public virtual void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (b.Length != orderA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (ReferenceEquals(a, b))
|
{
|
throw new ArgumentException("Arguments must be different objects.");
|
}
|
|
var clone = new Complex32[a.Length];
|
a.Copy(clone);
|
CholeskyFactor(clone, orderA);
|
CholeskySolveFactored(clone, orderA, b, columnsB);
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using a previously factored A matrix.
|
/// </summary>
|
/// <param name="a">The square, positive definite matrix A.</param>
|
/// <param name="orderA">The number of rows and columns in A.</param>
|
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
|
/// <param name="columnsB">The number of columns in the B matrix.</param>
|
/// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
|
public virtual void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (b.Length != orderA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (ReferenceEquals(a, b))
|
{
|
throw new ArgumentException("Arguments must be different objects.");
|
}
|
|
CommonParallel.For(0, columnsB, (u, v) =>
|
{
|
for (int i = u; i < v; i++)
|
{
|
DoCholeskySolve(a, orderA, b, i);
|
}
|
});
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using a previously factored A matrix.
|
/// </summary>
|
/// <param name="a">The square, positive definite matrix A. Has to be different than <paramref name="b"/>.</param>
|
/// <param name="orderA">The number of rows and columns in A.</param>
|
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
|
/// <param name="index">The column to solve for.</param>
|
static void DoCholeskySolve(Complex32[] a, int orderA, Complex32[] b, int index)
|
{
|
var cindex = index*orderA;
|
|
// Solve L*Y = B;
|
Complex32 sum;
|
for (var i = 0; i < orderA; i++)
|
{
|
sum = b[cindex + i];
|
for (var k = i - 1; k >= 0; k--)
|
{
|
sum -= a[(k*orderA) + i]*b[cindex + k];
|
}
|
|
b[cindex + i] = sum/a[(i*orderA) + i];
|
}
|
|
// Solve L'*X = Y;
|
for (var i = orderA - 1; i >= 0; i--)
|
{
|
sum = b[cindex + i];
|
var iindex = i*orderA;
|
for (var k = i + 1; k < orderA; k++)
|
{
|
sum -= a[iindex + k].Conjugate()*b[cindex + k];
|
}
|
|
b[cindex + i] = sum/a[iindex + i];
|
}
|
}
|
|
/// <summary>
|
/// Computes the QR factorization of A.
|
/// </summary>
|
/// <param name="r">On entry, it is the M by N A matrix to factor. On exit,
|
/// it is overwritten with the R matrix of the QR factorization. </param>
|
/// <param name="rowsR">The number of rows in the A matrix.</param>
|
/// <param name="columnsR">The number of columns in the A matrix.</param>
|
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
|
/// QR factorization.</param>
|
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
|
/// to be used by the QR solve routine.</param>
|
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
|
public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau)
|
{
|
if (r == null)
|
{
|
throw new ArgumentNullException(nameof(r));
|
}
|
|
if (q == null)
|
{
|
throw new ArgumentNullException(nameof(q));
|
}
|
|
if (r.Length != rowsR*columnsR)
|
{
|
throw new ArgumentException("The given array has the wrong length. Should be rowsR * columnsR.", nameof(r));
|
}
|
|
if (tau.Length < Math.Min(rowsR, columnsR))
|
{
|
throw new ArgumentException("The given array is too small. It must be at least min(m,n) long.", nameof(tau));
|
}
|
|
if (q.Length != rowsR*rowsR)
|
{
|
throw new ArgumentException("The given array has the wrong length. Should be rowsR * rowsR.", nameof(q));
|
}
|
|
var work = columnsR > rowsR ? new Complex32[rowsR*rowsR] : new Complex32[rowsR*columnsR];
|
|
CommonParallel.For(0, rowsR, (a, b) =>
|
{
|
for (int i = a; i < b; i++)
|
{
|
q[(i*rowsR) + i] = Complex32.One;
|
}
|
});
|
|
var minmn = Math.Min(rowsR, columnsR);
|
for (var i = 0; i < minmn; i++)
|
{
|
GenerateColumn(work, r, rowsR, i, i);
|
ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.MaxDegreeOfParallelism);
|
}
|
|
for (var i = minmn - 1; i >= 0; i--)
|
{
|
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.MaxDegreeOfParallelism);
|
}
|
}
|
|
/// <summary>
|
/// Computes the QR factorization of A.
|
/// </summary>
|
/// <param name="a">On entry, it is the M by N A matrix to factor. On exit,
|
/// it is overwritten with the Q matrix of the QR factorization.</param>
|
/// <param name="rowsA">The number of rows in the A matrix.</param>
|
/// <param name="columnsA">The number of columns in the A matrix.</param>
|
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
|
/// QR factorization.</param>
|
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
|
/// to be used by the QR solve routine.</param>
|
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
|
public virtual void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau)
|
{
|
if (r == null)
|
{
|
throw new ArgumentNullException(nameof(r));
|
}
|
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (a.Length != rowsA*columnsA)
|
{
|
throw new ArgumentException("The given array has the wrong length. Should be rowsR * columnsR.", nameof(a));
|
}
|
|
if (tau.Length < Math.Min(rowsA, columnsA))
|
{
|
throw new ArgumentException("The given array is too small. It must be at least min(m,n) long.", nameof(tau));
|
}
|
|
if (r.Length != columnsA*columnsA)
|
{
|
throw new ArgumentException("The given array has the wrong length. Should be columnsA * columnsA.", nameof(r));
|
}
|
|
var work = new Complex32[rowsA*columnsA];
|
|
var minmn = Math.Min(rowsA, columnsA);
|
for (var i = 0; i < minmn; i++)
|
{
|
GenerateColumn(work, a, rowsA, i, i);
|
ComputeQR(work, i, a, i, rowsA, i + 1, columnsA, Control.MaxDegreeOfParallelism);
|
}
|
|
//copy R
|
for (var j = 0; j < columnsA; j++)
|
{
|
var rIndex = j*columnsA;
|
var aIndex = j*rowsA;
|
for (var i = 0; i < columnsA; i++)
|
{
|
r[rIndex + i] = a[aIndex + i];
|
}
|
}
|
|
//clear A and set diagonals to 1
|
Array.Clear(a, 0, a.Length);
|
for (var i = 0; i < columnsA; i++)
|
{
|
a[i*rowsA + i] = Complex32.One;
|
}
|
|
for (var i = minmn - 1; i >= 0; i--)
|
{
|
ComputeQR(work, i, a, i, rowsA, i, columnsA, Control.MaxDegreeOfParallelism);
|
}
|
}
|
|
|
#region QR Factor Helper functions
|
|
/// <summary>
|
/// Perform calculation of Q or R
|
/// </summary>
|
/// <param name="work">Work array</param>
|
/// <param name="workIndex">Index of column in work array</param>
|
/// <param name="a">Q or R matrices</param>
|
/// <param name="rowStart">The first row in </param>
|
/// <param name="rowCount">The last row</param>
|
/// <param name="columnStart">The first column</param>
|
/// <param name="columnCount">The last column</param>
|
/// <param name="availableCores">Number of available CPUs</param>
|
static void ComputeQR(Complex32[] work, int workIndex, Complex32[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
|
{
|
if (rowStart > rowCount || columnStart > columnCount)
|
{
|
return;
|
}
|
|
var tmpColCount = columnCount - columnStart;
|
|
if ((availableCores > 1) && (tmpColCount > 200))
|
{
|
var tmpSplit = columnStart + (tmpColCount/2);
|
var tmpCores = availableCores/2;
|
|
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 < columnCount; j++)
|
{
|
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;
|
}
|
}
|
}
|
}
|
|
/// <summary>
|
/// Generate column from initial matrix to work array
|
/// </summary>
|
/// <param name="work">Work array</param>
|
/// <param name="a">Initial matrix</param>
|
/// <param name="rowCount">The number of rows in matrix</param>
|
/// <param name="row">The first row</param>
|
/// <param name="column">Column index</param>
|
static void GenerateColumn(Complex32[] work, Complex32[] a, int rowCount, int row, int column)
|
{
|
var tmp = column*rowCount;
|
var index = tmp + row;
|
|
CommonParallel.For(row, rowCount, (u, v) =>
|
{
|
for (int i = u; i < v; i++)
|
{
|
var iIndex = tmp + i;
|
work[iIndex - row] = a[iIndex];
|
a[iIndex] = Complex32.Zero;
|
}
|
});
|
|
var norm = Complex32.Zero;
|
for (var i = 0; i < rowCount - row; ++i)
|
{
|
var index1 = tmp + i;
|
norm += work[index1].Magnitude*work[index1].Magnitude;
|
}
|
|
norm = norm.SquareRoot();
|
if (row == rowCount - 1 || norm.Magnitude == 0)
|
{
|
a[index] = -work[tmp];
|
work[tmp] = new Complex32(2.0f, 0).SquareRoot();
|
return;
|
}
|
|
if (work[tmp].Magnitude != 0.0f)
|
{
|
norm = norm.Magnitude*(work[tmp]/work[tmp].Magnitude);
|
}
|
|
a[index] = -norm;
|
CommonParallel.For(0, rowCount - row, 4096, (u, v) =>
|
{
|
for (int i = u; i < v; i++)
|
{
|
work[tmp + i] /= norm;
|
}
|
});
|
work[tmp] += 1.0f;
|
|
var s = (1.0f/work[tmp]).SquareRoot();
|
CommonParallel.For(0, rowCount - row, 4096, (u, v) =>
|
{
|
for (int i = u; i < v; i++)
|
{
|
work[tmp + i] = work[tmp + i].Conjugate()*s;
|
}
|
});
|
}
|
|
#endregion
|
|
/// <summary>
|
/// Solves A*X=B for X using QR factorization of A.
|
/// </summary>
|
/// <param name="a">The A matrix.</param>
|
/// <param name="rows">The number of rows in the A matrix.</param>
|
/// <param name="columns">The number of columns in the A matrix.</param>
|
/// <param name="b">The B matrix.</param>
|
/// <param name="columnsB">The number of columns of B.</param>
|
/// <param name="x">On exit, the solution matrix.</param>
|
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
|
/// <remarks>Rows must be greater or equal to columns.</remarks>
|
public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (a.Length != rows*columns)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(a));
|
}
|
|
if (b.Length != rows*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (x.Length != columns*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(x));
|
}
|
|
if (rows < columns)
|
{
|
throw new ArgumentException("The number of rows must greater than or equal to the number of columns.");
|
}
|
|
var work = new Complex32[rows * columns];
|
|
var clone = new Complex32[a.Length];
|
a.Copy(clone);
|
|
if (method == QRMethod.Full)
|
{
|
var q = new Complex32[rows*rows];
|
QRFactor(clone, rows, columns, q, work);
|
QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x, method);
|
}
|
else
|
{
|
var r = new Complex32[columns*columns];
|
ThinQRFactor(clone, rows, columns, r, work);
|
QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
|
}
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using a previously QR factored matrix.
|
/// </summary>
|
/// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>.</param>
|
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
|
/// <param name="rowsA">The number of rows in the A matrix.</param>
|
/// <param name="columnsA">The number of columns in the A matrix.</param>
|
/// <param name="tau">Contains additional information on Q. Only used for the native solver
|
/// and can be <c>null</c> for the managed provider.</param>
|
/// <param name="b">The B matrix.</param>
|
/// <param name="columnsB">The number of columns of B.</param>
|
/// <param name="x">On exit, the solution matrix.</param>
|
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
|
/// <remarks>Rows must be greater or equal to columns.</remarks>
|
public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full)
|
{
|
if (r == null)
|
{
|
throw new ArgumentNullException(nameof(r));
|
}
|
|
if (q == null)
|
{
|
throw new ArgumentNullException(nameof(q));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(q));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(q));
|
}
|
|
if (rowsA < columnsA)
|
{
|
throw new ArgumentException("The number of rows must greater than or equal to the number of columns.");
|
}
|
|
int rowsQ, columnsQ, rowsR, columnsR;
|
if (method == QRMethod.Full)
|
{
|
rowsQ = columnsQ = rowsR = rowsA;
|
columnsR = columnsA;
|
}
|
else
|
{
|
rowsQ = rowsA;
|
columnsQ = rowsR = columnsR = columnsA;
|
}
|
|
if (r.Length != rowsR*columnsR)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {rowsR * columnsR}.", nameof(r));
|
}
|
|
if (q.Length != rowsQ*columnsQ)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {rowsQ * columnsQ}.", nameof(q));
|
}
|
|
if (b.Length != rowsA*columnsB)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {rowsA * columnsB}.", nameof(b));
|
}
|
|
if (x.Length != columnsA*columnsB)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {columnsA * columnsB}.", nameof(x));
|
}
|
|
var sol = new Complex32[b.Length];
|
|
// Copy B matrix to "sol", so B data will not be changed
|
Array.Copy(b, 0, sol, 0, b.Length);
|
|
// Compute Y = transpose(Q)*B
|
var column = new Complex32[rowsA];
|
for (var j = 0; j < columnsB; j++)
|
{
|
var jm = j*rowsA;
|
Array.Copy(sol, jm, column, 0, rowsA);
|
CommonParallel.For(0, columnsA, (u, v) =>
|
{
|
for (int i = u; i < v; i++)
|
{
|
var im = i*rowsA;
|
|
var sum = Complex32.Zero;
|
for (var k = 0; k < rowsA; k++)
|
{
|
sum += q[im + k].Conjugate()*column[k];
|
}
|
|
sol[jm + i] = sum;
|
}
|
});
|
}
|
|
// Solve R*X = Y;
|
for (var k = columnsA - 1; k >= 0; k--)
|
{
|
var km = k*rowsR;
|
for (var j = 0; j < columnsB; j++)
|
{
|
sol[(j*rowsA) + k] /= r[km + k];
|
}
|
|
for (var i = 0; i < k; i++)
|
{
|
for (var j = 0; j < columnsB; j++)
|
{
|
var jm = j*rowsA;
|
sol[jm + i] -= sol[jm + k]*r[km + i];
|
}
|
}
|
}
|
|
// Fill result matrix
|
for (var col = 0; col < columnsB; col++)
|
{
|
Array.Copy(sol, col*rowsA, x, col*columnsA, columnsR);
|
}
|
}
|
|
/// <summary>
|
/// Computes the singular value decomposition of A.
|
/// </summary>
|
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
|
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
|
/// <param name="rowsA">The number of rows in the A matrix.</param>
|
/// <param name="columnsA">The number of columns in the A matrix.</param>
|
/// <param name="s">The singular values of A in ascending value.</param>
|
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
|
/// singular vectors.</param>
|
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
|
/// right singular vectors.</param>
|
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
|
public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (s == null)
|
{
|
throw new ArgumentNullException(nameof(s));
|
}
|
|
if (u == null)
|
{
|
throw new ArgumentNullException(nameof(u));
|
}
|
|
if (vt == null)
|
{
|
throw new ArgumentNullException(nameof(vt));
|
}
|
|
if (u.Length != rowsA*rowsA)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(u));
|
}
|
|
if (vt.Length != columnsA*columnsA)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(vt));
|
}
|
|
if (s.Length != Math.Min(rowsA, columnsA))
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(s));
|
}
|
|
var work = new Complex32[rowsA];
|
|
const int maxiter = 1000;
|
|
var e = new Complex32[columnsA];
|
var v = new Complex32[vt.Length];
|
var stemp = new Complex32[Math.Min(rowsA + 1, columnsA)];
|
|
int i, j, l, lp1;
|
|
Complex32 t;
|
|
var ncu = rowsA;
|
|
// Reduce matrix to bidiagonal form, storing the diagonal elements
|
// in "s" and the super-diagonal elements in "e".
|
var nct = Math.Min(rowsA - 1, columnsA);
|
var nrt = Math.Max(0, Math.Min(columnsA - 2, rowsA));
|
var lu = Math.Max(nct, nrt);
|
|
for (l = 0; l < lu; l++)
|
{
|
lp1 = l + 1;
|
if (l < nct)
|
{
|
// Compute the transformation for the l-th column and
|
// place the l-th diagonal in vector s[l].
|
var sum = 0.0f;
|
for (i = l; i < rowsA; i++)
|
{
|
sum += a[(l*rowsA) + i].Magnitude*a[(l*rowsA) + i].Magnitude;
|
}
|
|
stemp[l] = (float) Math.Sqrt(sum);
|
if (stemp[l] != 0.0f)
|
{
|
if (a[(l*rowsA) + l] != 0.0f)
|
{
|
stemp[l] = stemp[l].Magnitude*(a[(l*rowsA) + l]/a[(l*rowsA) + l].Magnitude);
|
}
|
|
// A part of column "l" of Matrix A from row "l" to end multiply by 1.0f / s[l]
|
for (i = l; i < rowsA; i++)
|
{
|
a[(l*rowsA) + i] = a[(l*rowsA) + i]*(1.0f/stemp[l]);
|
}
|
|
a[(l*rowsA) + l] = 1.0f + a[(l*rowsA) + l];
|
}
|
|
stemp[l] = -stemp[l];
|
}
|
|
for (j = lp1; j < columnsA; j++)
|
{
|
if (l < nct)
|
{
|
if (stemp[l] != 0.0f)
|
{
|
// Apply the transformation.
|
t = 0.0f;
|
for (i = l; i < rowsA; i++)
|
{
|
t += a[(l*rowsA) + i].Conjugate()*a[(j*rowsA) + i];
|
}
|
|
t = -t/a[(l*rowsA) + l];
|
|
for (var ii = l; ii < rowsA; ii++)
|
{
|
a[(j*rowsA) + ii] += t*a[(l*rowsA) + ii];
|
}
|
}
|
}
|
|
// Place the l-th row of matrix into "e" for the
|
// subsequent calculation of the row transformation.
|
e[j] = a[(j*rowsA) + l].Conjugate();
|
}
|
|
if (computeVectors && l < nct)
|
{
|
// Place the transformation in "u" for subsequent back multiplication.
|
for (i = l; i < rowsA; i++)
|
{
|
u[(l*rowsA) + i] = a[(l*rowsA) + i];
|
}
|
}
|
|
if (l >= nrt)
|
{
|
continue;
|
}
|
|
// Compute the l-th row transformation and place the l-th super-diagonal in e(l).
|
var enorm = 0.0f;
|
for (i = lp1; i < e.Length; i++)
|
{
|
enorm += e[i].Magnitude*e[i].Magnitude;
|
}
|
|
e[l] = (float) Math.Sqrt(enorm);
|
if (e[l] != 0.0f)
|
{
|
if (e[lp1] != 0.0f)
|
{
|
e[l] = e[l].Magnitude*(e[lp1]/e[lp1].Magnitude);
|
}
|
|
// Scale vector "e" from "lp1" by 1.0f / e[l]
|
for (i = lp1; i < e.Length; i++)
|
{
|
e[i] = e[i]*(1.0f/e[l]);
|
}
|
|
e[lp1] = 1.0f + e[lp1];
|
}
|
|
e[l] = -e[l].Conjugate();
|
|
if (lp1 < rowsA && e[l] != 0.0f)
|
{
|
// Apply the transformation.
|
for (i = lp1; i < rowsA; i++)
|
{
|
work[i] = 0.0f;
|
}
|
|
for (j = lp1; j < columnsA; j++)
|
{
|
for (var ii = lp1; ii < rowsA; ii++)
|
{
|
work[ii] += e[j]*a[(j*rowsA) + ii];
|
}
|
}
|
|
for (j = lp1; j < columnsA; j++)
|
{
|
var ww = (-e[j]/e[lp1]).Conjugate();
|
for (var ii = lp1; ii < rowsA; ii++)
|
{
|
a[(j*rowsA) + ii] += ww*work[ii];
|
}
|
}
|
}
|
|
if (!computeVectors)
|
{
|
continue;
|
}
|
|
// Place the transformation in v for subsequent back multiplication.
|
for (i = lp1; i < columnsA; i++)
|
{
|
v[(l*columnsA) + i] = e[i];
|
}
|
}
|
|
// Set up the final bidiagonal matrix or order m.
|
var m = Math.Min(columnsA, rowsA + 1);
|
var nctp1 = nct + 1;
|
var nrtp1 = nrt + 1;
|
if (nct < columnsA)
|
{
|
stemp[nctp1 - 1] = a[((nctp1 - 1)*rowsA) + (nctp1 - 1)];
|
}
|
|
if (rowsA < m)
|
{
|
stemp[m - 1] = 0.0f;
|
}
|
|
if (nrtp1 < m)
|
{
|
e[nrtp1 - 1] = a[((m - 1)*rowsA) + (nrtp1 - 1)];
|
}
|
|
e[m - 1] = 0.0f;
|
|
// If required, generate "u".
|
if (computeVectors)
|
{
|
for (j = nctp1 - 1; j < ncu; j++)
|
{
|
for (i = 0; i < rowsA; i++)
|
{
|
u[(j*rowsA) + i] = 0.0f;
|
}
|
|
u[(j*rowsA) + j] = 1.0f;
|
}
|
|
for (l = nct - 1; l >= 0; l--)
|
{
|
if (stemp[l] != 0.0f)
|
{
|
for (j = l + 1; j < ncu; j++)
|
{
|
t = 0.0f;
|
for (i = l; i < rowsA; i++)
|
{
|
t += u[(l*rowsA) + i].Conjugate()*u[(j*rowsA) + i];
|
}
|
|
t = -t/u[(l*rowsA) + l];
|
for (var ii = l; ii < rowsA; ii++)
|
{
|
u[(j*rowsA) + ii] += t*u[(l*rowsA) + ii];
|
}
|
}
|
|
// A part of column "l" of matrix A from row "l" to end multiply by -1.0f
|
for (i = l; i < rowsA; i++)
|
{
|
u[(l*rowsA) + i] = u[(l*rowsA) + i]*-1.0f;
|
}
|
|
u[(l*rowsA) + l] = 1.0f + u[(l*rowsA) + l];
|
for (i = 0; i < l; i++)
|
{
|
u[(l*rowsA) + i] = 0.0f;
|
}
|
}
|
else
|
{
|
for (i = 0; i < rowsA; i++)
|
{
|
u[(l*rowsA) + i] = 0.0f;
|
}
|
|
u[(l*rowsA) + l] = 1.0f;
|
}
|
}
|
}
|
|
// If it is required, generate v.
|
if (computeVectors)
|
{
|
for (l = columnsA - 1; l >= 0; l--)
|
{
|
lp1 = l + 1;
|
if (l < nrt)
|
{
|
if (e[l] != 0.0f)
|
{
|
for (j = lp1; j < columnsA; j++)
|
{
|
t = 0.0f;
|
for (i = lp1; i < columnsA; i++)
|
{
|
t += v[(l*columnsA) + i].Conjugate()*v[(j*columnsA) + i];
|
}
|
|
t = -t/v[(l*columnsA) + lp1];
|
for (var ii = l; ii < columnsA; ii++)
|
{
|
v[(j*columnsA) + ii] += t*v[(l*columnsA) + ii];
|
}
|
}
|
}
|
}
|
|
for (i = 0; i < columnsA; i++)
|
{
|
v[(l*columnsA) + i] = 0.0f;
|
}
|
|
v[(l*columnsA) + l] = 1.0f;
|
}
|
}
|
|
// Transform "s" and "e" so that they are float
|
for (i = 0; i < m; i++)
|
{
|
Complex32 r;
|
if (stemp[i] != 0.0f)
|
{
|
t = stemp[i].Magnitude;
|
r = stemp[i]/t;
|
stemp[i] = t;
|
if (i < m - 1)
|
{
|
e[i] = e[i]/r;
|
}
|
|
if (computeVectors)
|
{
|
// A part of column "i" of matrix U from row 0 to end multiply by r
|
for (j = 0; j < rowsA; j++)
|
{
|
u[(i*rowsA) + j] = u[(i*rowsA) + j]*r;
|
}
|
}
|
}
|
|
// Exit
|
if (i == m - 1)
|
{
|
break;
|
}
|
|
if (e[i] == 0.0f)
|
{
|
continue;
|
}
|
|
t = e[i].Magnitude;
|
r = t/e[i];
|
e[i] = t;
|
stemp[i + 1] = stemp[i + 1]*r;
|
if (!computeVectors)
|
{
|
continue;
|
}
|
|
// A part of column "i+1" of matrix VT from row 0 to end multiply by r
|
for (j = 0; j < columnsA; j++)
|
{
|
v[((i + 1)*columnsA) + j] = v[((i + 1)*columnsA) + j]*r;
|
}
|
}
|
|
// Main iteration loop for the singular values.
|
var mn = m;
|
var iter = 0;
|
|
while (m > 0)
|
{
|
// Quit if all the singular values have been found.
|
// If too many iterations have been performed throw exception.
|
if (iter >= maxiter)
|
{
|
throw new NonConvergenceException();
|
}
|
|
// This section of the program inspects for negligible elements in the s and e arrays,
|
// on completion the variables case and l are set as follows:
|
// case = 1: if mS[m] and e[l-1] are negligible and l < m
|
// case = 2: if mS[l] is negligible and l < m
|
// case = 3: if e[l-1] is negligible, l < m, and mS[l, ..., mS[m] are not negligible (qr step).
|
// case = 4: if e[m-1] is negligible (convergence).
|
float ztest;
|
float test;
|
for (l = m - 2; l >= 0; l--)
|
{
|
test = stemp[l].Magnitude + stemp[l + 1].Magnitude;
|
ztest = test + e[l].Magnitude;
|
if (ztest.AlmostEqualRelative(test, 7))
|
{
|
e[l] = 0.0f;
|
break;
|
}
|
}
|
|
int kase;
|
if (l == m - 2)
|
{
|
kase = 4;
|
}
|
else
|
{
|
int ls;
|
for (ls = m - 1; ls > l; ls--)
|
{
|
test = 0.0f;
|
if (ls != m - 1)
|
{
|
test = test + e[ls].Magnitude;
|
}
|
|
if (ls != l + 1)
|
{
|
test = test + e[ls - 1].Magnitude;
|
}
|
|
ztest = test + stemp[ls].Magnitude;
|
if (ztest.AlmostEqualRelative(test, 7))
|
{
|
stemp[ls] = 0.0f;
|
break;
|
}
|
}
|
|
if (ls == l)
|
{
|
kase = 3;
|
}
|
else if (ls == m - 1)
|
{
|
kase = 1;
|
}
|
else
|
{
|
kase = 2;
|
l = ls;
|
}
|
}
|
|
l = l + 1;
|
|
// Perform the task indicated by case.
|
int k;
|
float f;
|
float sn;
|
float cs;
|
switch (kase)
|
{
|
// Deflate negligible s[m].
|
case 1:
|
f = e[m - 2].Real;
|
e[m - 2] = 0.0f;
|
float t1;
|
for (var kk = l; kk < m - 1; kk++)
|
{
|
k = m - 2 - kk + l;
|
t1 = stemp[k].Real;
|
Drotg(ref t1, ref f, out cs, out sn);
|
stemp[k] = t1;
|
if (k != l)
|
{
|
f = -sn*e[k - 1].Real;
|
e[k - 1] = cs*e[k - 1];
|
}
|
|
if (computeVectors)
|
{
|
// Rotate
|
for (i = 0; i < columnsA; i++)
|
{
|
var z = (cs*v[(k*columnsA) + i]) + (sn*v[((m - 1)*columnsA) + i]);
|
v[((m - 1)*columnsA) + i] = (cs*v[((m - 1)*columnsA) + i]) - (sn*v[(k*columnsA) + i]);
|
v[(k*columnsA) + i] = z;
|
}
|
}
|
}
|
|
break;
|
|
// Split at negligible s[l].
|
case 2:
|
f = e[l - 1].Real;
|
e[l - 1] = 0.0f;
|
for (k = l; k < m; k++)
|
{
|
t1 = stemp[k].Real;
|
Drotg(ref t1, ref f, out cs, out sn);
|
stemp[k] = t1;
|
f = -sn*e[k].Real;
|
e[k] = cs*e[k];
|
if (computeVectors)
|
{
|
// Rotate
|
for (i = 0; i < rowsA; i++)
|
{
|
var z = (cs*u[(k*rowsA) + i]) + (sn*u[((l - 1)*rowsA) + i]);
|
u[((l - 1)*rowsA) + i] = (cs*u[((l - 1)*rowsA) + i]) - (sn*u[(k*rowsA) + i]);
|
u[(k*rowsA) + i] = z;
|
}
|
}
|
}
|
|
break;
|
|
// Perform one qr step.
|
case 3:
|
// calculate the shift.
|
var scale = 0.0f;
|
scale = Math.Max(scale, stemp[m - 1].Magnitude);
|
scale = Math.Max(scale, stemp[m - 2].Magnitude);
|
scale = Math.Max(scale, e[m - 2].Magnitude);
|
scale = Math.Max(scale, stemp[l].Magnitude);
|
scale = Math.Max(scale, e[l].Magnitude);
|
var sm = stemp[m - 1].Real/scale;
|
var smm1 = stemp[m - 2].Real/scale;
|
var emm1 = e[m - 2].Real/scale;
|
var sl = stemp[l].Real/scale;
|
var el = e[l].Real/scale;
|
var b = (((smm1 + sm)*(smm1 - sm)) + (emm1*emm1))/2.0f;
|
var c = (sm*emm1)*(sm*emm1);
|
var shift = 0.0f;
|
if (b != 0.0f || c != 0.0f)
|
{
|
shift = (float) Math.Sqrt((b*b) + c);
|
if (b < 0.0f)
|
{
|
shift = -shift;
|
}
|
|
shift = c/(b + shift);
|
}
|
|
f = ((sl + sm)*(sl - sm)) + shift;
|
var g = sl*el;
|
|
// Chase zeros
|
for (k = l; k < m - 1; k++)
|
{
|
Drotg(ref f, ref g, out cs, out sn);
|
if (k != l)
|
{
|
e[k - 1] = f;
|
}
|
|
f = (cs*stemp[k].Real) + (sn*e[k].Real);
|
e[k] = (cs*e[k]) - (sn*stemp[k]);
|
g = sn*stemp[k + 1].Real;
|
stemp[k + 1] = cs*stemp[k + 1];
|
if (computeVectors)
|
{
|
for (i = 0; i < columnsA; i++)
|
{
|
var z = (cs*v[(k*columnsA) + i]) + (sn*v[((k + 1)*columnsA) + i]);
|
v[((k + 1)*columnsA) + i] = (cs*v[((k + 1)*columnsA) + i]) - (sn*v[(k*columnsA) + i]);
|
v[(k*columnsA) + i] = z;
|
}
|
}
|
|
Drotg(ref f, ref g, out cs, out sn);
|
stemp[k] = f;
|
f = (cs*e[k].Real) + (sn*stemp[k + 1].Real);
|
stemp[k + 1] = -(sn*e[k]) + (cs*stemp[k + 1]);
|
g = sn*e[k + 1].Real;
|
e[k + 1] = cs*e[k + 1];
|
if (computeVectors && k < rowsA)
|
{
|
for (i = 0; i < rowsA; i++)
|
{
|
var z = (cs*u[(k*rowsA) + i]) + (sn*u[((k + 1)*rowsA) + i]);
|
u[((k + 1)*rowsA) + i] = (cs*u[((k + 1)*rowsA) + i]) - (sn*u[(k*rowsA) + i]);
|
u[(k*rowsA) + i] = z;
|
}
|
}
|
}
|
|
e[m - 2] = f;
|
iter = iter + 1;
|
break;
|
|
// Convergence
|
case 4:
|
|
// Make the singular value positive
|
if (stemp[l].Real < 0.0f)
|
{
|
stemp[l] = -stemp[l];
|
if (computeVectors)
|
{
|
// A part of column "l" of matrix VT from row 0 to end multiply by -1
|
for (i = 0; i < columnsA; i++)
|
{
|
v[(l*columnsA) + i] = v[(l*columnsA) + i]*-1.0f;
|
}
|
}
|
}
|
|
// Order the singular value.
|
while (l != mn - 1)
|
{
|
if (stemp[l].Real >= stemp[l + 1].Real)
|
{
|
break;
|
}
|
|
t = stemp[l];
|
stemp[l] = stemp[l + 1];
|
stemp[l + 1] = t;
|
if (computeVectors && l < columnsA)
|
{
|
// Swap columns l, l + 1
|
for (i = 0; i < columnsA; i++)
|
{
|
var z = v[(l*columnsA) + i];
|
v[(l*columnsA) + i] = v[((l + 1)*columnsA) + i];
|
v[((l + 1)*columnsA) + i] = z;
|
}
|
}
|
|
if (computeVectors && l < rowsA)
|
{
|
// Swap columns l, l + 1
|
for (i = 0; i < rowsA; i++)
|
{
|
var z = u[(l*rowsA) + i];
|
u[(l*rowsA) + i] = u[((l + 1)*rowsA) + i];
|
u[((l + 1)*rowsA) + i] = z;
|
}
|
}
|
|
l = l + 1;
|
}
|
|
iter = 0;
|
m = m - 1;
|
break;
|
}
|
}
|
|
if (computeVectors)
|
{
|
// Finally transpose "v" to get "vt" matrix
|
for (i = 0; i < columnsA; i++)
|
{
|
for (j = 0; j < columnsA; j++)
|
{
|
vt[(j*columnsA) + i] = v[(i*columnsA) + j].Conjugate();
|
}
|
}
|
}
|
|
// Copy stemp to s with size adjustment. We are using ported copy of linpack's svd code and it uses
|
// a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed.
|
// We should port lapack's svd routine to remove this problem.
|
Array.Copy(stemp, 0, s, 0, Math.Min(rowsA, columnsA));
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using the singular value decomposition of A.
|
/// </summary>
|
/// <param name="a">On entry, the M by N matrix to decompose.</param>
|
/// <param name="rowsA">The number of rows in the A matrix.</param>
|
/// <param name="columnsA">The number of columns in the A matrix.</param>
|
/// <param name="b">The B matrix.</param>
|
/// <param name="columnsB">The number of columns of B.</param>
|
/// <param name="x">On exit, the solution matrix.</param>
|
public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x)
|
{
|
if (a == null)
|
{
|
throw new ArgumentNullException(nameof(a));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (b.Length != rowsA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (x.Length != columnsA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
var s = new Complex32[Math.Min(rowsA, columnsA)];
|
var u = new Complex32[rowsA*rowsA];
|
var vt = new Complex32[columnsA*columnsA];
|
|
var clone = new Complex32[a.Length];
|
a.Copy(clone);
|
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
|
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
|
}
|
|
/// <summary>
|
/// Solves A*X=B for X using a previously SVD decomposed matrix.
|
/// </summary>
|
/// <param name="rowsA">The number of rows in the A matrix.</param>
|
/// <param name="columnsA">The number of columns in the A matrix.</param>
|
/// <param name="s">The s values returned by <see cref="SingularValueDecomposition(bool,Complex32[],int,int,Complex32[],Complex32[],Complex32[])"/>.</param>
|
/// <param name="u">The left singular vectors returned by <see cref="SingularValueDecomposition(bool,Complex32[],int,int,Complex32[],Complex32[],Complex32[])"/>.</param>
|
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,Complex32[],int,int,Complex32[],Complex32[],Complex32[])"/>.</param>
|
/// <param name="b">The B matrix.</param>
|
/// <param name="columnsB">The number of columns of B.</param>
|
/// <param name="x">On exit, the solution matrix.</param>
|
public virtual void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x)
|
{
|
if (s == null)
|
{
|
throw new ArgumentNullException(nameof(s));
|
}
|
|
if (u == null)
|
{
|
throw new ArgumentNullException(nameof(u));
|
}
|
|
if (vt == null)
|
{
|
throw new ArgumentNullException(nameof(vt));
|
}
|
|
if (b == null)
|
{
|
throw new ArgumentNullException(nameof(b));
|
}
|
|
if (x == null)
|
{
|
throw new ArgumentNullException(nameof(x));
|
}
|
|
if (u.Length != rowsA*rowsA)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(u));
|
}
|
|
if (vt.Length != columnsA*columnsA)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(vt));
|
}
|
|
if (s.Length != Math.Min(rowsA, columnsA))
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(s));
|
}
|
|
if (b.Length != rowsA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
if (x.Length != columnsA*columnsB)
|
{
|
throw new ArgumentException("The array arguments must have the same length.", nameof(b));
|
}
|
|
var mn = Math.Min(rowsA, columnsA);
|
var tmp = new Complex32[columnsA];
|
|
for (var k = 0; k < columnsB; k++)
|
{
|
for (var j = 0; j < columnsA; j++)
|
{
|
var value = Complex32.Zero;
|
if (j < mn)
|
{
|
for (var i = 0; i < rowsA; i++)
|
{
|
value += u[(j*rowsA) + i].Conjugate()*b[(k*rowsA) + i];
|
}
|
|
value /= s[j];
|
}
|
|
tmp[j] = value;
|
}
|
|
for (var j = 0; j < columnsA; j++)
|
{
|
var value = Complex32.Zero;
|
for (var i = 0; i < columnsA; i++)
|
{
|
value += vt[(j*columnsA) + i].Conjugate()*tmp[i];
|
}
|
|
x[(k*columnsA) + j] = value;
|
}
|
}
|
}
|
|
/// <summary>
|
/// Computes the eigenvalues and eigenvectors of a matrix.
|
/// </summary>
|
/// <param name="isSymmetric">Whether the matrix is symmetric or not.</param>
|
/// <param name="order">The order of the matrix.</param>
|
/// <param name="matrix">The matrix to decompose. The length of the array must be order * order.</param>
|
/// <param name="matrixEv">On output, the matrix contains the eigen vectors. The length of the array must be order * order.</param>
|
/// <param name="vectorEv">On output, the eigen values (λ) of matrix in ascending value. The length of the array must <paramref name="order"/>.</param>
|
/// <param name="matrixD">On output, the block diagonal eigenvalue matrix. The length of the array must be order * order.</param>
|
public virtual void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD)
|
{
|
if (matrix == null)
|
{
|
throw new ArgumentNullException(nameof(matrix));
|
}
|
|
if (matrix.Length != order*order)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {order * order}.", nameof(matrix));
|
}
|
|
if (matrixEv == null)
|
{
|
throw new ArgumentNullException(nameof(matrixEv));
|
}
|
|
if (matrixEv.Length != order*order)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {order * order}.", nameof(matrixEv));
|
}
|
|
if (vectorEv == null)
|
{
|
throw new ArgumentNullException(nameof(vectorEv));
|
}
|
|
if (vectorEv.Length != order)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {order}.", nameof(vectorEv));
|
}
|
|
if (matrixD == null)
|
{
|
throw new ArgumentNullException(nameof(matrixD));
|
}
|
|
if (matrixD.Length != order*order)
|
{
|
throw new ArgumentException($"The given array has the wrong length. Should be {order * order}.", nameof(matrixD));
|
}
|
|
var matrixCopy = new Complex32[matrix.Length];
|
Array.Copy(matrix, 0, matrixCopy, 0, matrix.Length);
|
if (isSymmetric)
|
{
|
var tau = new Complex32[order];
|
var d = new float[order];
|
var e = new float[order];
|
|
Managed.ManagedLinearAlgebraProvider.SymmetricTridiagonalize(matrixCopy, d, e, tau, order);
|
Managed.ManagedLinearAlgebraProvider.SymmetricDiagonalize(matrixEv, d, e, order);
|
Managed.ManagedLinearAlgebraProvider.SymmetricUntridiagonalize(matrixEv, matrixCopy, tau, order);
|
|
for (var i = 0; i < order; i++)
|
{
|
vectorEv[i] = new Complex(d[i], e[i]);
|
matrixD[i*order + i] = new Complex32(d[i], e[i]);
|
}
|
}
|
else
|
{
|
var v = new Complex32[order];
|
|
Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceToHessenberg(matrixEv, matrixCopy, order);
|
Managed.ManagedLinearAlgebraProvider.NonsymmetricReduceHessenberToRealSchur(v, matrixEv, matrixCopy, order);
|
|
for (var i = 0; i < order; i++)
|
{
|
vectorEv[i] = new Complex(v[i].Real, v[i].Imaginary);
|
matrixD[i*order + i] = v[i];
|
}
|
}
|
}
|
|
/// <summary>
|
/// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
|
/// </summary>
|
protected static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, Complex32[] matrix, Complex32[] row)
|
{
|
if (transpose == Transpose.DontTranspose)
|
{
|
for (int i = 0; i < numCols; i++)
|
{
|
row[i] = matrix[(i*numRows) + rowindx];
|
}
|
}
|
else if (transpose == Transpose.ConjugateTranspose)
|
{
|
int offset = rowindx*numCols;
|
for (int i = 0; i < row.Length; i++)
|
{
|
row[i] = matrix[i + offset].Conjugate();
|
}
|
}
|
else
|
{
|
Array.Copy(matrix, rowindx*numCols, row, 0, numCols);
|
}
|
}
|
|
/// <summary>
|
/// Assumes that <paramref name="numRows"/> and <paramref name="numCols"/> have already been transposed.
|
/// </summary>
|
protected static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, Complex32[] matrix, Complex32[] column)
|
{
|
if (transpose == Transpose.DontTranspose)
|
{
|
Array.Copy(matrix, colindx*numRows, column, 0, numRows);
|
}
|
else if (transpose == Transpose.ConjugateTranspose)
|
{
|
for (int i = 0; i < numRows; i++)
|
{
|
column[i] = matrix[(i*numCols) + colindx].Conjugate();
|
}
|
}
|
else
|
{
|
for (int i = 0; i < numRows; i++)
|
{
|
column[i] = matrix[(i*numCols) + colindx];
|
}
|
}
|
}
|
}
|
}
|
