//
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
// Copyright (c) 2009-2010 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.
//
using System.Numerics;
using MathNet.Numerics.LinearAlgebra.Generic.Factorization;
namespace MathNet.Numerics.Algorithms.LinearAlgebra
{
using System;
using Properties;
using Threading;
///
/// The managed linear algebra provider.
///
public partial class ManagedLinearAlgebraProvider
{
///
/// Adds a scaled vector to another: result = y + alpha*x.
///
/// The vector to update.
/// The value to scale by.
/// The vector to add to .
/// The result of the addition.
/// This is similar to the AXPY BLAS routine.
public virtual void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (y.Length != x.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (alpha.IsZero())
{
y.Copy(result);
}
else if (alpha.IsOne())
{
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = y[i] + x[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = y[index] + x[index];
}
}
}
else
{
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = y[i] + (alpha*x[i]);
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = y[index] + (alpha * x[index]);
}
}
}
}
///
/// Scales an array. Can be used to scale a vector and a matrix.
///
/// The scalar.
/// The values to scale.
/// This result of the scaling.
/// This is similar to the SCAL BLAS routine.
public virtual void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result)
{
if (x == null)
{
throw new ArgumentNullException("x");
}
if (alpha.IsZero())
{
Array.Clear(result, 0, result.Length);
}
else if (alpha.IsOne())
{
x.Copy(result);
}
else
{
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, x.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = alpha*x[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = alpha * x[index];
}
}
}
}
///
/// Computes the dot product of x and y.
///
/// The vector x.
/// The vector y.
/// The dot product of x and y.
/// This is equivalent to the DOT BLAS routine.
public virtual Complex32 DotProduct(Complex32[] x, Complex32[] y)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (y.Length != x.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
var d = new Complex32(0.0F, 0.0F);
for (var i = 0; i < y.Length; i++)
{
d += y[i] * x[i];
}
return d;
}
///
/// Does a point wise add of two arrays z = x + y. This can be used
/// to add vectors or matrices.
///
/// The array x.
/// The array y.
/// The result of the addition.
/// There is no equivalent BLAS routine, but many libraries
/// provide optimized (parallel and/or vectorized) versions of this
/// routine.
public virtual void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (y.Length != x.Length || y.Length != result.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = x[i] + y[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = x[index] + y[index];
}
}
}
///
/// Does a point wise subtraction of two arrays z = x - y. This can be used
/// to subtract vectors or matrices.
///
/// The array x.
/// The array y.
/// The result of the subtraction.
/// There is no equivalent BLAS routine, but many libraries
/// provide optimized (parallel and/or vectorized) versions of this
/// routine.
public virtual void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (y.Length != x.Length || y.Length != result.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = x[i] - y[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = x[index] - y[index];
}
}
}
///
/// Does a point wise multiplication of two arrays z = x * y. This can be used
/// to multiple elements of vectors or matrices.
///
/// The array x.
/// The array y.
/// The result of the point wise multiplication.
/// There is no equivalent BLAS routine, but many libraries
/// provide optimized (parallel and/or vectorized) versions of this
/// routine.
public virtual void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (y.Length != x.Length || y.Length != result.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = x[i]*y[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = x[index] * y[index];
}
}
}
///
/// Does a point wise division of two arrays z = x / y. This can be used
/// to divide elements of vectors or matrices.
///
/// The array x.
/// The array y.
/// The result of the point wise division.
/// There is no equivalent BLAS routine, but many libraries
/// provide optimized (parallel and/or vectorized) versions of this
/// routine.
public virtual void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result)
{
if (y == null)
{
throw new ArgumentNullException("y");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (result == null)
{
throw new ArgumentNullException("result");
}
if (y.Length != x.Length || y.Length != result.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (Control.ParallelizeOperation(x.Length))
{
CommonParallel.For(0, y.Length, (a, b) =>
{
for (int i = a; i < b; i++)
{
result[i] = x[i]/y[i];
}
});
}
else
{
for (var index = 0; index < x.Length; index++)
{
result[index] = x[index] / y[index];
}
}
}
///
/// Computes the requested of the matrix.
///
/// The type of norm to compute.
/// The number of rows.
/// The number of columns.
/// The matrix to compute the norm from.
///
/// The requested of the matrix.
///
public virtual Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix)
{
var ret = 0.0;
switch (norm)
{
case Norm.OneNorm:
for (var j = 0; j < columns; j++)
{
var s = 0.0;
for (var i = 0; i < rows; i++)
{
s += matrix[(j * rows) + i].Magnitude;
}
ret = Math.Max(ret, s);
}
break;
case Norm.LargestAbsoluteValue:
for (var i = 0; i < rows; i++)
{
for (var j = 0; j < columns; j++)
{
ret = Math.Max(matrix[(j * rows) + i].Magnitude, ret);
}
}
break;
case Norm.InfinityNorm:
for (var i = 0; i < rows; i++)
{
var s = 0.0;
for (var j = 0; j < columns; j++)
{
s += matrix[(j * rows) + i].Magnitude;
}
ret = Math.Max(ret, s);
}
break;
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);
for (var i = 0; i < rows; i++)
{
ret += aat[(i * rows) + i].Magnitude;
}
ret = Math.Sqrt(ret);
break;
}
return Convert.ToSingle(ret);
}
///
/// Computes the requested of the matrix.
///
/// The type of norm to compute.
/// The number of rows.
/// The number of columns.
/// The matrix to compute the norm from.
/// The work array. Only used when
/// and needs to be have a length of at least M (number of rows of .
///
/// The requested of the matrix.
///
public virtual Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix, float[] work)
{
return MatrixNorm(norm, rows, columns, matrix);
}
///
/// Multiples two matrices. result = x * y
///
/// The x matrix.
/// The number of rows in the x matrix.
/// The number of columns in the x matrix.
/// The y matrix.
/// The number of rows in the y matrix.
/// The number of columns in the y matrix.
/// Where to store the result of the multiplication.
/// 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.
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("x");
}
if (y == null)
{
throw new ArgumentNullException("y");
}
if (result == null)
{
throw new ArgumentNullException("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;
}
MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, xdata, rowsX, columnsX, ydata, rowsY, columnsY, Complex32.Zero, result);
}
///
/// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c
///
/// How to transpose the matrix.
/// How to transpose the matrix.
/// The value to scale matrix.
/// The a matrix.
/// The number of rows in the matrix.
/// The number of columns in the matrix.
/// The b matrix
/// The number of rows in the matrix.
/// The number of columns in the matrix.
/// The value to scale the matrix.
/// The c matrix.
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("a");
}
if (b == null)
{
throw new ArgumentNullException("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())
{
Control.LinearAlgebraProvider.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);
}
///
/// Cache-Oblivious Matrix Multiplication
///
/// if set to true transpose matrix A.
/// if set to true transpose matrix B.
/// The value to scale the matrix A with.
/// The matrix A.
/// Row-shift of the left matrix
/// Column-shift of the left matrix
/// The matrix B.
/// Row-shift of the right matrix
/// Column-shift of the right matrix
/// The matrix C.
/// Row-shift of the result matrix
/// Column-shift of the result matrix
/// The number of rows of matrix op(A) and of the matrix C.
/// The number of columns of matrix op(B) and of the matrix C.
/// The number of columns of matrix op(A) and the rows of the matrix op(B).
/// The constant number of rows of matrix op(A) and of the matrix C.
/// The constant number of columns of matrix op(B) and of the matrix C.
/// The constant number of columns of matrix op(A) and the rows of the matrix op(B).
/// Indicates if this is the first recursion.
private 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)
{
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);
}
}
}
///
/// Computes the LUP factorization of A. P*A = L*U.
///
/// An by matrix. The matrix is overwritten with the
/// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0
/// for the L factor). The upper triangular factor U is stored on and above the diagonal of .
/// The order of the square matrix .
/// On exit, it contains the pivot indices. The size of the array must be .
/// This is equivalent to the GETRF LAPACK routine.
public virtual void LUFactor(Complex32[] data, int order, int[] ipiv)
{
if (data == null)
{
throw new ArgumentNullException("data");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (data.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "data");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "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];
}
}
}
}
///
/// Computes the inverse of matrix using LU factorization.
///
/// The N by N matrix to invert. Contains the inverse On exit.
/// The order of the square matrix .
/// This is equivalent to the GETRF and GETRI LAPACK routines.
public virtual void LUInverse(Complex32[] a, int order)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
var ipiv = new int[order];
LUFactor(a, order, ipiv);
LUInverseFactored(a, order, ipiv);
}
///
/// Computes the inverse of a previously factored matrix.
///
/// The LU factored N by N matrix. Contains the inverse On exit.
/// The order of the square matrix .
/// The pivot indices of .
/// This is equivalent to the GETRI LAPACK routine.
public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "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);
}
///
/// Computes the inverse of matrix using LU factorization.
///
/// The N by N matrix to invert. Contains the inverse On exit.
/// The order of the square matrix .
/// The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// This is equivalent to the GETRF and GETRI LAPACK routines.
public virtual void LUInverse(Complex32[] a, int order, Complex32[] work)
{
LUInverse(a, order);
}
///
/// Computes the inverse of a previously factored matrix.
///
/// The LU factored N by N matrix. Contains the inverse On exit.
/// The order of the square matrix .
/// The pivot indices of .
/// The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// This is equivalent to the GETRI LAPACK routine.
public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work)
{
LUInverseFactored(a, order, ipiv);
}
///
/// Solves A*X=B for X using LU factorization.
///
/// The number of columns of B.
/// The square matrix A.
/// The order of the square matrix .
/// On entry the B matrix; on exit the X matrix.
/// This is equivalent to the GETRF and GETRS LAPACK routines.
public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (ReferenceEquals(a, b))
{
throw new ArgumentException(Resources.ArgumentReferenceDifferent);
}
var ipiv = new int[order];
var clone = new Complex32[a.Length];
a.Copy(clone);
LUFactor(clone, order, ipiv);
LUSolveFactored(columnsOfB, clone, order, ipiv, b);
}
///
/// Solves A*X=B for X using a previously factored A matrix.
///
/// The number of columns of B.
/// The factored A matrix.
/// The order of the square matrix .
/// The pivot indices of .
/// On entry the B matrix; on exit the X matrix.
/// This is equivalent to the GETRS LAPACK routine.
public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (a.Length != order * order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (b.Length != order * columnsOfB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (ReferenceEquals(a, b))
{
throw new ArgumentException(Resources.ArgumentReferenceDifferent);
}
// 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];
}
}
}
}
///
/// Computes the Cholesky factorization of A.
///
/// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the
/// the Cholesky factorization.
/// The number of rows or columns in the matrix.
/// This is equivalent to the POTRF LAPACK routine.
public virtual void CholeskyFactor(Complex32[] a, int order)
{
if (a == null)
{
throw new ArgumentNullException("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.NumberOfParallelWorkerThreads);
}
else
{
throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
}
for (var i = ij + 1; i < order; i++)
{
a[(i * order) + ij] = 0.0f;
}
}
}
///
/// Calculate Cholesky step
///
/// Factor matrix
/// Number of rows
/// Column start
/// Total columns
/// Multipliers calculated previously
/// Number of available processors
private 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();
}
}
}
}
///
/// Solves A*X=B for X using Cholesky factorization.
///
/// The square, positive definite matrix A.
/// The number of rows and columns in A.
/// On entry the B matrix; on exit the X matrix.
/// The number of columns in the B matrix.
/// This is equivalent to the POTRF add POTRS LAPACK routines.
public virtual void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (b.Length != orderA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (ReferenceEquals(a, b))
{
throw new ArgumentException(Resources.ArgumentReferenceDifferent);
}
var clone = new Complex32[a.Length];
a.Copy(clone);
CholeskyFactor(clone, orderA);
CholeskySolveFactored(clone, orderA, b, columnsB);
}
///
/// Solves A*X=B for X using a previously factored A matrix.
///
/// The square, positive definite matrix A.
/// The number of rows and columns in A.
/// On entry the B matrix; on exit the X matrix.
/// The number of columns in the B matrix.
/// This is equivalent to the POTRS LAPACK routine.
public virtual void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (b.Length != orderA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (ReferenceEquals(a, b))
{
throw new ArgumentException(Resources.ArgumentReferenceDifferent);
}
CommonParallel.For(0, columnsB, (u, v) =>
{
for (int i = u; i < v; i++)
{
DoCholeskySolve(a, orderA, b, i);
}
});
}
///
/// Solves A*X=B for X using a previously factored A matrix.
///
/// The square, positive definite matrix A. Has to be different than .
/// The number of rows and columns in A.
/// On entry the B matrix; on exit the X matrix.
/// The column to solve for.
private 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];
}
}
///
/// Computes the QR factorization of A.
///
/// 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.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.
/// A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.
/// This is similar to the GEQRF and ORGQR LAPACK routines.
public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR * columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR * rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
var work = columnsR > rowsR ? new Complex32 [rowsR * rowsR] : new Complex32[rowsR * columnsR];
QRFactor(r, rowsR, columnsR, q, tau, work);
}
///
/// Computes the QR factorization of A.
///
/// 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.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.
/// A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.
/// The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// This is similar to the GEQRF and ORGQR LAPACK routines.
public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau, Complex32[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR * columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR * rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
if (columnsR > rowsR)
{
if (work.Length < rowsR * rowsR)
{
work[0] = rowsR * rowsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
else
{
if (work.Length < rowsR * columnsR)
{
work[0] = rowsR * columnsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
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.NumberOfParallelWorkerThreads);
}
for (var i = minmn - 1; i >= 0; i--)
{
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
}
work[0] = columnsR > rowsR ? rowsR * rowsR : rowsR * columnsR;
}
///
/// Computes the QR factorization of 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.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.
/// A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.
/// This is similar to the GEQRF and ORGQR LAPACK routines.
public virtual void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != rowsA * columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA * columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
var work = new Complex32[rowsA * columnsA];
ThinQRFactor(a, rowsA, columnsA, r, tau, work);
}
///
/// Computes the QR factorization of A where M > N.
///
/// 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.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.
/// A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.
/// The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// This is similar to the GEQRF and ORGQR LAPACK routines.
public virtual void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau, Complex32[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (a.Length != rowsA * columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA * columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
if (work.Length < rowsA * columnsA)
{
work[0] = rowsA * columnsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
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.NumberOfParallelWorkerThreads);
}
//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.NumberOfParallelWorkerThreads);
}
work[0] = rowsA * columnsA;
}
#region QR Factor Helper functions
///
/// Perform calculation of Q or R
///
/// Work array
/// Index of column in work array
/// Q or R matrices
/// The first row in
/// The last row
/// The first column
/// The last column
/// Number of available CPUs
private 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;
}
}
}
}
///
/// Generate column from initial matrix to work array
///
/// Work array
/// Initial matrix
/// The number of rows in matrix
/// The first row
/// Column index
private 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, (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, (u, v) =>
{
for (int i = u; i < v; i++)
{
work[tmp + i] = work[tmp + i].Conjugate()*s;
}
});
}
#endregion
///
/// Solves A*X=B for X using QR factorization of A.
///
/// The A matrix.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
/// The type of QR factorization to perform.
/// Rows must be greater or equal to columns.
public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full)
{
var work = new Complex32[rows * columns];
QRSolve(a, rows, columns, b, columnsB, x, work, method);
}
///
/// Solves A*X=B for X using QR factorization of A.
///
/// The A matrix.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
/// The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// The type of QR factorization to perform.
/// Rows must be greater or equal to columns.
public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rows * columns)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (b.Length != rows * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (x.Length != columns * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
}
if (rows < columns)
{
throw new ArgumentException(Resources.RowsLessThanColumns);
}
if (work.Length < rows * columns)
{
work[0] = rows * columns;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
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);
}
work[0] = rows * columns;
}
///
/// Solves A*X=B for X using a previously QR factored matrix.
///
/// The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// null for the native provider. The native provider uses the Q portion stored in the R matrix.
/// The R matrix obtained by calling .
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// Contains additional information on Q. Only used for the native solver
/// and can be null for the managed provider.
/// On entry the B matrix; on exit the X matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
/// The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.
/// The type of QR factorization to perform.
/// Rows must be greater or equal to columns.
public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full)
{
QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, method);
}
///
/// Solves A*X=B for X using a previously QR factored matrix.
///
/// The Q matrix obtained by calling .
/// The R matrix obtained by calling .
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// Contains additional information on Q. Only used for the native solver
/// and can be null for the managed provider.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
/// The type of QR factorization to perform.
/// Rows must be greater or equal to columns.
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("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (b == null)
{
throw new ArgumentNullException("q");
}
if (x == null)
{
throw new ArgumentNullException("q");
}
if (rowsA < columnsA)
{
throw new ArgumentException(Resources.RowsLessThanColumns);
}
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(string.Format(Resources.ArgumentArrayWrongLength, rowsR*columnsR), "r");
}
if (q.Length != rowsQ*columnsQ)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ*columnsQ), "q");
}
if (b.Length != rowsA*columnsB)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA*columnsB), "b");
}
if (x.Length != columnsA*columnsB)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA*columnsB), "x");
}
var sol = new Complex32[b.Length];
// Copy B matrix to "sol", so B data will not be changed
Array.Copy(b, sol, b.Length);
// Compute Y = transpose(Q)*B
var column = new Complex32[rowsA];
for (var j = 0; j < columnsB; j++)
{
var jm = j*rowsA;
CommonParallel.For(0, rowsA, (u, v) =>
{
for (int k = u; k < v; k++)
{
column[k] = sol[jm + k];
}
});
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
CommonParallel.For(0, columnsR, (u, v) =>
{
for (int row = u; row < v; row++)
{
for (var col = 0; col < columnsB; col++)
{
x[(col*columnsA) + row] = sol[row + (col*rowsA)];
}
}
});
}
///
/// Computes the singular value decomposition of A.
///
/// Compute the singular U and VT vectors or not.
/// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The singular values of A in ascending value.
/// If is true, on exit U contains the left
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (u.Length != rowsA * rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA * columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
var work = new Complex32[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
///
/// Computes the singular value decomposition of A.
///
/// Compute the singular U and VT vectors or not.
/// On entry, the M by N matrix to decompose. On exit, A may be overwritten.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The singular values of A in ascending value.
/// If is true, on exit U contains the left
/// singular vectors.
/// If is true, on exit VT contains the transposed
/// right singular vectors.
/// The work array. Length should be at least .
/// This is equivalent to the GESVD LAPACK routine.
public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (u.Length != rowsA * rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA * columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (work.Length == 0)
{
throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
}
if (work.Length < rowsA)
{
work[0] = rowsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
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;
var cs = 0.0f;
var sn = 0.0f;
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 ArgumentException(Resources.ConvergenceFailed);
}
// This section of the program inspects for negligible elements in the s and e arrays,
// on completion the variables kase and l are set as follows:
// kase = 1: if mS[m] and e[l-1] are negligible and l < m
// kase = 2: if mS[l] is negligible and l < m
// kase = 3: if e[l-1] is negligible, l < m, and mS[l, ..., mS[m] are not negligible (qr step).
// kase = 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.AlmostEqualInDecimalPlaces(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.AlmostEqualInDecimalPlaces(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 kase.
int k;
float f;
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, ref cs, ref 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, ref cs, ref 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, ref cs, ref 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, ref cs, ref 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, s, Math.Min(rowsA, columnsA));
// On return the first element of the work array stores the min size of the work array could have been
// work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns));
work[0] = rowsA;
}
///
/// Solves A*X=B for X using the singular value decomposition of A.
///
/// On entry, the M by N matrix to decompose.
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (x.Length != columnsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var work = new Complex32[rowsA];
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, work);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
///
/// Solves A*X=B for X using a previously SVD decomposed matrix.
///
/// The number of rows in the A matrix.
/// The number of columns in the A matrix.
/// The s values returned by .
/// The left singular vectors returned by .
/// The right singular vectors returned by .
/// The B matrix.
/// The number of columns of B.
/// On exit, the solution matrix.
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("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (b == null)
{
throw new ArgumentNullException("b");
}
if (x == null)
{
throw new ArgumentNullException("x");
}
if (u.Length != rowsA * rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA * columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (b.Length != rowsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
if (x.Length != columnsA * columnsB)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "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;
}
}
}
}
}