|
|
|
@ -4,7 +4,7 @@ |
|
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
|
//
|
|
|
|
// Copyright (c) 2009-2010 Math.NET
|
|
|
|
// Copyright (c) 2009-2013 Math.NET
|
|
|
|
//
|
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
|
@ -28,8 +28,8 @@ |
|
|
|
// OTHER DEALINGS IN THE SOFTWARE.
|
|
|
|
// </copyright>
|
|
|
|
|
|
|
|
using MathNet.Numerics.Properties; |
|
|
|
using System; |
|
|
|
using MathNet.Numerics.Properties; |
|
|
|
|
|
|
|
namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
{ |
|
|
|
@ -49,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <remarks>
|
|
|
|
/// The computation of the singular value decomposition is done at construction time.
|
|
|
|
/// </remarks>
|
|
|
|
public class UserSvd : Svd |
|
|
|
public sealed class UserSvd : Svd |
|
|
|
{ |
|
|
|
/// <summary>
|
|
|
|
/// Initializes a new instance of the <see cref="UserSvd"/> class. This object will compute the
|
|
|
|
@ -59,20 +59,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
|
|
|
|
/// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
|
|
|
|
/// <exception cref="NonConvergenceException"></exception>
|
|
|
|
public UserSvd(Matrix<Complex32> matrix, bool computeVectors) |
|
|
|
public static UserSvd Create(Matrix<Complex32> matrix, bool computeVectors) |
|
|
|
{ |
|
|
|
if (matrix == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("matrix"); |
|
|
|
} |
|
|
|
|
|
|
|
ComputeVectors = computeVectors; |
|
|
|
var nm = Math.Min(matrix.RowCount + 1, matrix.ColumnCount); |
|
|
|
var matrixCopy = matrix.Clone(); |
|
|
|
|
|
|
|
S = matrixCopy.CreateVector(nm); |
|
|
|
U = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); |
|
|
|
VT = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); |
|
|
|
var s = matrixCopy.CreateVector(nm); |
|
|
|
var u = matrixCopy.CreateMatrix(matrixCopy.RowCount, matrixCopy.RowCount); |
|
|
|
var vt = matrixCopy.CreateMatrix(matrixCopy.ColumnCount, matrixCopy.ColumnCount); |
|
|
|
|
|
|
|
const int maxiter = 1000; |
|
|
|
var e = new Complex32[matrixCopy.ColumnCount]; |
|
|
|
@ -95,34 +89,34 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
if (l < nct) |
|
|
|
{ |
|
|
|
// Compute the transformation for the l-th column and place the l-th diagonal in VectorS[l].
|
|
|
|
S[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); |
|
|
|
if (S[l].Magnitude != 0.0f) |
|
|
|
s[l] = Cnrm2Column(matrixCopy, matrixCopy.RowCount, l, l); |
|
|
|
if (s[l].Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
if (matrixCopy.At(l, l).Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
S[l] = Csign(S[l], matrixCopy.At(l, l)); |
|
|
|
s[l] = Csign(s[l], matrixCopy.At(l, l)); |
|
|
|
} |
|
|
|
|
|
|
|
CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f / S[l]); |
|
|
|
CscalColumn(matrixCopy, matrixCopy.RowCount, l, l, 1.0f/s[l]); |
|
|
|
matrixCopy.At(l, l, (Complex32.One + matrixCopy.At(l, l))); |
|
|
|
} |
|
|
|
|
|
|
|
S[l] = -S[l]; |
|
|
|
s[l] = -s[l]; |
|
|
|
} |
|
|
|
|
|
|
|
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|
|
|
{ |
|
|
|
if (l < nct) |
|
|
|
{ |
|
|
|
if (S[l].Magnitude != 0.0f) |
|
|
|
if (s[l].Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
// Apply the transformation.
|
|
|
|
t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l) / matrixCopy.At(l, l); |
|
|
|
t = -Cdotc(matrixCopy, matrixCopy.RowCount, l, j, l)/matrixCopy.At(l, l); |
|
|
|
if (t != Complex32.Zero) |
|
|
|
{ |
|
|
|
for (var ii = l; ii < matrixCopy.RowCount; ii++) |
|
|
|
{ |
|
|
|
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t * matrixCopy.At(ii, l))); |
|
|
|
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (t*matrixCopy.At(ii, l))); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
@ -133,12 +127,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
e[j] = matrixCopy.At(l, j).Conjugate(); |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors && l < nct) |
|
|
|
if (computeVectors && l < nct) |
|
|
|
{ |
|
|
|
// Place the transformation in u for subsequent back multiplication.
|
|
|
|
for (i = l; i < matrixCopy.RowCount; i++) |
|
|
|
{ |
|
|
|
U.At(i, l, matrixCopy.At(i, l)); |
|
|
|
u.At(i, l, matrixCopy.At(i, l)); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -157,7 +151,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
e[l] = Csign(e[l], e[lp1]); |
|
|
|
} |
|
|
|
|
|
|
|
CscalVector(e, lp1, 1.0f / e[l]); |
|
|
|
CscalVector(e, lp1, 1.0f/e[l]); |
|
|
|
e[lp1] = Complex32.One + e[lp1]; |
|
|
|
} |
|
|
|
|
|
|
|
@ -176,30 +170,30 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
{ |
|
|
|
for (var ii = lp1; ii < matrixCopy.RowCount; ii++) |
|
|
|
{ |
|
|
|
work[ii] += e[j] * matrixCopy.At(ii, j); |
|
|
|
work[ii] += e[j]*matrixCopy.At(ii, j); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|
|
|
{ |
|
|
|
var ww = (-e[j] / e[lp1]).Conjugate(); |
|
|
|
var ww = (-e[j]/e[lp1]).Conjugate(); |
|
|
|
if (ww != Complex32.Zero) |
|
|
|
{ |
|
|
|
for (var ii = lp1; ii < matrixCopy.RowCount; ii++) |
|
|
|
{ |
|
|
|
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww * work[ii])); |
|
|
|
matrixCopy.At(ii, j, matrixCopy.At(ii, j) + (ww*work[ii])); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
// Place the transformation in v for subsequent back multiplication.
|
|
|
|
for (i = lp1; i < matrixCopy.ColumnCount; i++) |
|
|
|
{ |
|
|
|
VT.At(i, l, e[i]); |
|
|
|
vt.At(i, l, e[i]); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
@ -210,12 +204,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
var nrtp1 = nrt + 1; |
|
|
|
if (nct < matrixCopy.ColumnCount) |
|
|
|
{ |
|
|
|
S[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); |
|
|
|
s[nctp1 - 1] = matrixCopy.At((nctp1 - 1), (nctp1 - 1)); |
|
|
|
} |
|
|
|
|
|
|
|
if (matrixCopy.RowCount < m) |
|
|
|
{ |
|
|
|
S[m - 1] = Complex32.Zero; |
|
|
|
s[m - 1] = Complex32.Zero; |
|
|
|
} |
|
|
|
|
|
|
|
if (nrtp1 < m) |
|
|
|
@ -226,55 +220,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
e[m - 1] = Complex32.Zero; |
|
|
|
|
|
|
|
// If required, generate u.
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
for (j = nctp1 - 1; j < ncu; j++) |
|
|
|
{ |
|
|
|
for (i = 0; i < matrixCopy.RowCount; i++) |
|
|
|
{ |
|
|
|
U.At(i, j, Complex32.Zero); |
|
|
|
u.At(i, j, Complex32.Zero); |
|
|
|
} |
|
|
|
|
|
|
|
U.At(j, j, Complex32.One); |
|
|
|
u.At(j, j, Complex32.One); |
|
|
|
} |
|
|
|
|
|
|
|
for (l = nct - 1; l >= 0; l--) |
|
|
|
{ |
|
|
|
if (S[l].Magnitude != 0.0f) |
|
|
|
if (s[l].Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
for (j = l + 1; j < ncu; j++) |
|
|
|
{ |
|
|
|
t = -Cdotc(U, matrixCopy.RowCount, l, j, l) / U.At(l, l); |
|
|
|
t = -Cdotc(u, matrixCopy.RowCount, l, j, l)/u.At(l, l); |
|
|
|
if (t != Complex32.Zero) |
|
|
|
{ |
|
|
|
for (var ii = l; ii < matrixCopy.RowCount; ii++) |
|
|
|
{ |
|
|
|
U.At(ii, j, U.At(ii, j) + (t * U.At(ii, l))); |
|
|
|
u.At(ii, j, u.At(ii, j) + (t*u.At(ii, l))); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
CscalColumn(U, matrixCopy.RowCount, l, l, -1.0f); |
|
|
|
U.At(l, l, Complex32.One + U.At(l, l)); |
|
|
|
CscalColumn(u, matrixCopy.RowCount, l, l, -1.0f); |
|
|
|
u.At(l, l, Complex32.One + u.At(l, l)); |
|
|
|
for (i = 0; i < l; i++) |
|
|
|
{ |
|
|
|
U.At(i, l, Complex32.Zero); |
|
|
|
u.At(i, l, Complex32.Zero); |
|
|
|
} |
|
|
|
} |
|
|
|
else |
|
|
|
{ |
|
|
|
for (i = 0; i < matrixCopy.RowCount; i++) |
|
|
|
{ |
|
|
|
U.At(i, l, Complex32.Zero); |
|
|
|
u.At(i, l, Complex32.Zero); |
|
|
|
} |
|
|
|
|
|
|
|
U.At(l, l, Complex32.One); |
|
|
|
u.At(l, l, Complex32.One); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
// If it is required, generate v.
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
for (l = matrixCopy.ColumnCount - 1; l >= 0; l--) |
|
|
|
{ |
|
|
|
@ -285,12 +279,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
{ |
|
|
|
for (j = lp1; j < matrixCopy.ColumnCount; j++) |
|
|
|
{ |
|
|
|
t = -Cdotc(VT, matrixCopy.ColumnCount, l, j, lp1) / VT.At(lp1, l); |
|
|
|
t = -Cdotc(vt, matrixCopy.ColumnCount, l, j, lp1)/vt.At(lp1, l); |
|
|
|
if (t != Complex32.Zero) |
|
|
|
{ |
|
|
|
for (var ii = l; ii < matrixCopy.ColumnCount; ii++) |
|
|
|
{ |
|
|
|
VT.At(ii, j, VT.At(ii, j) + (t * VT.At(ii, l))); |
|
|
|
vt.At(ii, j, vt.At(ii, j) + (t*vt.At(ii, l))); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
@ -299,10 +293,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
|
|
|
|
for (i = 0; i < matrixCopy.ColumnCount; i++) |
|
|
|
{ |
|
|
|
VT.At(i, l, Complex32.Zero); |
|
|
|
vt.At(i, l, Complex32.Zero); |
|
|
|
} |
|
|
|
|
|
|
|
VT.At(l, l, Complex32.One); |
|
|
|
vt.At(l, l, Complex32.One); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -310,19 +304,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
for (i = 0; i < m; i++) |
|
|
|
{ |
|
|
|
Complex32 r; |
|
|
|
if (S[i].Magnitude != 0.0f) |
|
|
|
if (s[i].Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
t = S[i].Magnitude; |
|
|
|
r = S[i] / t; |
|
|
|
S[i] = t; |
|
|
|
t = s[i].Magnitude; |
|
|
|
r = s[i]/t; |
|
|
|
s[i] = t; |
|
|
|
if (i < m - 1) |
|
|
|
{ |
|
|
|
e[i] = e[i] / r; |
|
|
|
e[i] = e[i]/r; |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
CscalColumn(U, matrixCopy.RowCount, i, 0, r); |
|
|
|
CscalColumn(u, matrixCopy.RowCount, i, 0, r); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -335,12 +329,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
if (e[i].Magnitude != 0.0f) |
|
|
|
{ |
|
|
|
t = e[i].Magnitude; |
|
|
|
r = t / e[i]; |
|
|
|
r = t/e[i]; |
|
|
|
e[i] = t; |
|
|
|
S[i + 1] = S[i + 1] * r; |
|
|
|
if (ComputeVectors) |
|
|
|
s[i + 1] = s[i + 1]*r; |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
CscalColumn(VT, matrixCopy.ColumnCount, i + 1, 0, r); |
|
|
|
CscalColumn(vt, matrixCopy.ColumnCount, i + 1, 0, r); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
@ -368,7 +362,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
float test; |
|
|
|
for (l = m - 2; l >= 0; l--) |
|
|
|
{ |
|
|
|
test = S[l].Magnitude + S[l + 1].Magnitude; |
|
|
|
test = s[l].Magnitude + s[l + 1].Magnitude; |
|
|
|
ztest = test + e[l].Magnitude; |
|
|
|
if (ztest.AlmostEqualInDecimalPlaces(test, 7)) |
|
|
|
{ |
|
|
|
@ -398,10 +392,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
test = test + e[ls - 1].Magnitude; |
|
|
|
} |
|
|
|
|
|
|
|
ztest = test + S[ls].Magnitude; |
|
|
|
ztest = test + s[ls].Magnitude; |
|
|
|
if (ztest.AlmostEqualInDecimalPlaces(test, 7)) |
|
|
|
{ |
|
|
|
S[ls] = Complex32.Zero; |
|
|
|
s[ls] = Complex32.Zero; |
|
|
|
break; |
|
|
|
} |
|
|
|
} |
|
|
|
@ -430,7 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
float cs; |
|
|
|
switch (kase) |
|
|
|
{ |
|
|
|
// Deflate negligible VectorS[m].
|
|
|
|
// Deflate negligible VectorS[m].
|
|
|
|
case 1: |
|
|
|
f = e[m - 2].Real; |
|
|
|
e[m - 2] = Complex32.Zero; |
|
|
|
@ -438,73 +432,73 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
for (var kk = l; kk < m - 1; kk++) |
|
|
|
{ |
|
|
|
k = m - 2 - kk + l; |
|
|
|
t1 = S[k].Real; |
|
|
|
t1 = s[k].Real; |
|
|
|
Srotg(ref t1, ref f, out cs, out sn); |
|
|
|
S[k] = t1; |
|
|
|
s[k] = t1; |
|
|
|
if (k != l) |
|
|
|
{ |
|
|
|
f = -sn * e[k - 1].Real; |
|
|
|
e[k - 1] = cs * e[k - 1]; |
|
|
|
f = -sn*e[k - 1].Real; |
|
|
|
e[k - 1] = cs*e[k - 1]; |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
Csrot(VT, matrixCopy.ColumnCount, k, m - 1, cs, sn); |
|
|
|
Csrot(vt, matrixCopy.ColumnCount, k, m - 1, cs, sn); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
break; |
|
|
|
|
|
|
|
// Split at negligible VectorS[l].
|
|
|
|
// Split at negligible VectorS[l].
|
|
|
|
case 2: |
|
|
|
f = e[l - 1].Real; |
|
|
|
e[l - 1] = Complex32.Zero; |
|
|
|
for (k = l; k < m; k++) |
|
|
|
{ |
|
|
|
t1 = S[k].Real; |
|
|
|
t1 = s[k].Real; |
|
|
|
Srotg(ref t1, ref f, out cs, out sn); |
|
|
|
S[k] = t1; |
|
|
|
f = -sn * e[k].Real; |
|
|
|
e[k] = cs * e[k]; |
|
|
|
if (ComputeVectors) |
|
|
|
s[k] = t1; |
|
|
|
f = -sn*e[k].Real; |
|
|
|
e[k] = cs*e[k]; |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
Csrot(U, matrixCopy.RowCount, k, l - 1, cs, sn); |
|
|
|
Csrot(u, matrixCopy.RowCount, k, l - 1, cs, sn); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
break; |
|
|
|
|
|
|
|
// Perform one qr step.
|
|
|
|
// Perform one qr step.
|
|
|
|
case 3: |
|
|
|
// Calculate the shift.
|
|
|
|
var scale = 0.0f; |
|
|
|
scale = Math.Max(scale, S[m - 1].Magnitude); |
|
|
|
scale = Math.Max(scale, S[m - 2].Magnitude); |
|
|
|
scale = Math.Max(scale, s[m - 1].Magnitude); |
|
|
|
scale = Math.Max(scale, s[m - 2].Magnitude); |
|
|
|
scale = Math.Max(scale, e[m - 2].Magnitude); |
|
|
|
scale = Math.Max(scale, S[l].Magnitude); |
|
|
|
scale = Math.Max(scale, s[l].Magnitude); |
|
|
|
scale = Math.Max(scale, e[l].Magnitude); |
|
|
|
var sm = S[m - 1].Real / scale; |
|
|
|
var smm1 = S[m - 2].Real / scale; |
|
|
|
var emm1 = e[m - 2].Real / scale; |
|
|
|
var sl = S[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 sm = s[m - 1].Real/scale; |
|
|
|
var smm1 = s[m - 2].Real/scale; |
|
|
|
var emm1 = e[m - 2].Real/scale; |
|
|
|
var sl = s[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); |
|
|
|
shift = (float) Math.Sqrt((b*b) + c); |
|
|
|
if (b < 0.0f) |
|
|
|
{ |
|
|
|
shift = -shift; |
|
|
|
} |
|
|
|
|
|
|
|
shift = c / (b + shift); |
|
|
|
shift = c/(b + shift); |
|
|
|
} |
|
|
|
|
|
|
|
f = ((sl + sm) * (sl - sm)) + shift; |
|
|
|
var g = sl * el; |
|
|
|
f = ((sl + sm)*(sl - sm)) + shift; |
|
|
|
var g = sl*el; |
|
|
|
|
|
|
|
// Chase zeros.
|
|
|
|
for (k = l; k < m - 1; k++) |
|
|
|
@ -515,24 +509,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
e[k - 1] = f; |
|
|
|
} |
|
|
|
|
|
|
|
f = (cs * S[k].Real) + (sn * e[k].Real); |
|
|
|
e[k] = (cs * e[k]) - (sn * S[k]); |
|
|
|
g = sn * S[k + 1].Real; |
|
|
|
S[k + 1] = cs * S[k + 1]; |
|
|
|
if (ComputeVectors) |
|
|
|
f = (cs*s[k].Real) + (sn*e[k].Real); |
|
|
|
e[k] = (cs*e[k]) - (sn*s[k]); |
|
|
|
g = sn*s[k + 1].Real; |
|
|
|
s[k + 1] = cs*s[k + 1]; |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
Csrot(VT, matrixCopy.ColumnCount, k, k + 1, cs, sn); |
|
|
|
Csrot(vt, matrixCopy.ColumnCount, k, k + 1, cs, sn); |
|
|
|
} |
|
|
|
|
|
|
|
Srotg(ref f, ref g, out cs, out sn); |
|
|
|
S[k] = f; |
|
|
|
f = (cs * e[k].Real) + (sn * S[k + 1].Real); |
|
|
|
S[k + 1] = (-sn * e[k]) + (cs * S[k + 1]); |
|
|
|
g = sn * e[k + 1].Real; |
|
|
|
e[k + 1] = cs * e[k + 1]; |
|
|
|
if (ComputeVectors && k < matrixCopy.RowCount) |
|
|
|
s[k] = f; |
|
|
|
f = (cs*e[k].Real) + (sn*s[k + 1].Real); |
|
|
|
s[k + 1] = (-sn*e[k]) + (cs*s[k + 1]); |
|
|
|
g = sn*e[k + 1].Real; |
|
|
|
e[k + 1] = cs*e[k + 1]; |
|
|
|
if (computeVectors && k < matrixCopy.RowCount) |
|
|
|
{ |
|
|
|
Csrot(U, matrixCopy.RowCount, k, k + 1, cs, sn); |
|
|
|
Csrot(u, matrixCopy.RowCount, k, k + 1, cs, sn); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -540,37 +534,37 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
iter = iter + 1; |
|
|
|
break; |
|
|
|
|
|
|
|
// Convergence.
|
|
|
|
// Convergence.
|
|
|
|
case 4: |
|
|
|
// Make the singular value positive
|
|
|
|
if (S[l].Real < 0.0f) |
|
|
|
if (s[l].Real < 0.0f) |
|
|
|
{ |
|
|
|
S[l] = -S[l]; |
|
|
|
if (ComputeVectors) |
|
|
|
s[l] = -s[l]; |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
CscalColumn(VT, matrixCopy.ColumnCount, l, 0, -1.0f); |
|
|
|
CscalColumn(vt, matrixCopy.ColumnCount, l, 0, -1.0f); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
// Order the singular value.
|
|
|
|
while (l != mn - 1) |
|
|
|
{ |
|
|
|
if (S[l].Real >= S[l + 1].Real) |
|
|
|
if (s[l].Real >= s[l + 1].Real) |
|
|
|
{ |
|
|
|
break; |
|
|
|
} |
|
|
|
|
|
|
|
t = S[l]; |
|
|
|
S[l] = S[l + 1]; |
|
|
|
S[l + 1] = t; |
|
|
|
if (ComputeVectors && l < matrixCopy.ColumnCount) |
|
|
|
t = s[l]; |
|
|
|
s[l] = s[l + 1]; |
|
|
|
s[l + 1] = t; |
|
|
|
if (computeVectors && l < matrixCopy.ColumnCount) |
|
|
|
{ |
|
|
|
Swap(VT, matrixCopy.ColumnCount, l, l + 1); |
|
|
|
Swap(vt, matrixCopy.ColumnCount, l, l + 1); |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors && l < matrixCopy.RowCount) |
|
|
|
if (computeVectors && l < matrixCopy.RowCount) |
|
|
|
{ |
|
|
|
Swap(U, matrixCopy.RowCount, l, l + 1); |
|
|
|
Swap(u, matrixCopy.RowCount, l, l + 1); |
|
|
|
} |
|
|
|
|
|
|
|
l = l + 1; |
|
|
|
@ -582,9 +576,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
if (ComputeVectors) |
|
|
|
if (computeVectors) |
|
|
|
{ |
|
|
|
VT = VT.ConjugateTranspose(); |
|
|
|
vt = vt.ConjugateTranspose(); |
|
|
|
} |
|
|
|
|
|
|
|
// Adjust the size of s if rows < columns. We are using ported copy of linpack's svd code and it uses
|
|
|
|
@ -596,11 +590,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
var tmp = matrixCopy.CreateVector(nm); |
|
|
|
for (i = 0; i < nm; i++) |
|
|
|
{ |
|
|
|
tmp[i] = S[i]; |
|
|
|
tmp[i] = s[i]; |
|
|
|
} |
|
|
|
|
|
|
|
S = tmp; |
|
|
|
s = tmp; |
|
|
|
} |
|
|
|
|
|
|
|
return new UserSvd(s, u, vt, computeVectors); |
|
|
|
} |
|
|
|
|
|
|
|
UserSvd(Vector<Complex32> s, Matrix<Complex32> u, Matrix<Complex32> vt, bool vectorsComputed) |
|
|
|
: base(s, u, vt, vectorsComputed) |
|
|
|
{ |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
@ -609,9 +610,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="z1">Complex32 value z1</param>
|
|
|
|
/// <param name="z2">Complex32 value z2</param>
|
|
|
|
/// <returns>Result multiplication of signum function and absolute value</returns>
|
|
|
|
private static Complex32 Csign(Complex32 z1, Complex32 z2) |
|
|
|
static Complex32 Csign(Complex32 z1, Complex32 z2) |
|
|
|
{ |
|
|
|
return z1.Magnitude * (z2 / z2.Magnitude); |
|
|
|
return z1.Magnitude*(z2/z2.Magnitude); |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
@ -621,7 +622,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="rowCount">The number of rows in <paramref name="a"/></param>
|
|
|
|
/// <param name="columnA">Column A index to swap</param>
|
|
|
|
/// <param name="columnB">Column B index to swap</param>
|
|
|
|
private static void Swap(Matrix<Complex32> a, int rowCount, int columnA, int columnB) |
|
|
|
static void Swap(Matrix<Complex32> a, int rowCount, int columnA, int columnB) |
|
|
|
{ |
|
|
|
for (var i = 0; i < rowCount; i++) |
|
|
|
{ |
|
|
|
@ -639,11 +640,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="column">Column to scale</param>
|
|
|
|
/// <param name="rowStart">Row to scale from</param>
|
|
|
|
/// <param name="z">Scale value</param>
|
|
|
|
private static void CscalColumn(Matrix<Complex32> a, int rowCount, int column, int rowStart, Complex32 z) |
|
|
|
static void CscalColumn(Matrix<Complex32> a, int rowCount, int column, int rowStart, Complex32 z) |
|
|
|
{ |
|
|
|
for (var i = rowStart; i < rowCount; i++) |
|
|
|
{ |
|
|
|
a.At(i, column, a.At(i, column) * z); |
|
|
|
a.At(i, column, a.At(i, column)*z); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -653,11 +654,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="a">Source vector</param>
|
|
|
|
/// <param name="start">Row to scale from</param>
|
|
|
|
/// <param name="z">Scale value</param>
|
|
|
|
private static void CscalVector(Complex32[] a, int start, Complex32 z) |
|
|
|
static void CscalVector(Complex32[] a, int start, Complex32 z) |
|
|
|
{ |
|
|
|
for (var i = start; i < a.Length; i++) |
|
|
|
{ |
|
|
|
a[i] = a[i] * z; |
|
|
|
a[i] = a[i]*z; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -670,7 +671,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="c">Contains the parameter c associated with the Givens rotation</param>
|
|
|
|
/// <param name="s">Contains the parameter s associated with the Givens rotation</param>
|
|
|
|
/// <remarks>This is equivalent to the DROTG LAPACK routine.</remarks>
|
|
|
|
private static void Srotg(ref float da, ref float db, out float c, out float s) |
|
|
|
static void Srotg(ref float da, ref float db, out float c, out float s) |
|
|
|
{ |
|
|
|
float r, z; |
|
|
|
|
|
|
|
@ -692,16 +693,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
} |
|
|
|
else |
|
|
|
{ |
|
|
|
var sda = da / scale; |
|
|
|
var sdb = db / scale; |
|
|
|
r = scale * (float)Math.Sqrt((sda * sda) + (sdb * sdb)); |
|
|
|
var sda = da/scale; |
|
|
|
var sdb = db/scale; |
|
|
|
r = scale*(float) Math.Sqrt((sda*sda) + (sdb*sdb)); |
|
|
|
if (roe < 0.0f) |
|
|
|
{ |
|
|
|
r = -r; |
|
|
|
} |
|
|
|
|
|
|
|
c = da / r; |
|
|
|
s = db / r; |
|
|
|
c = da/r; |
|
|
|
s = db/r; |
|
|
|
z = 1.0f; |
|
|
|
if (absda > absdb) |
|
|
|
{ |
|
|
|
@ -710,7 +711,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
|
|
|
|
if (absdb >= absda && c != 0.0f) |
|
|
|
{ |
|
|
|
z = 1.0f / c; |
|
|
|
z = 1.0f/c; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
@ -726,15 +727,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="column">Column index</param>
|
|
|
|
/// <param name="rowStart">Start row index</param>
|
|
|
|
/// <returns>Norm2 (Euclidean norm) of the column</returns>
|
|
|
|
private static float Cnrm2Column(Matrix<Complex32> a, int rowCount, int column, int rowStart) |
|
|
|
static float Cnrm2Column(Matrix<Complex32> a, int rowCount, int column, int rowStart) |
|
|
|
{ |
|
|
|
var s = 0.0f; |
|
|
|
for (var i = rowStart; i < rowCount; i++) |
|
|
|
{ |
|
|
|
s += a.At(i, column).Magnitude * a.At(i, column).Magnitude; |
|
|
|
s += a.At(i, column).Magnitude*a.At(i, column).Magnitude; |
|
|
|
} |
|
|
|
|
|
|
|
return (float)Math.Sqrt(s); |
|
|
|
return (float) Math.Sqrt(s); |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
@ -743,15 +744,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="a">Source vector</param>
|
|
|
|
/// <param name="rowStart">Start index</param>
|
|
|
|
/// <returns>Norm2 (Euclidean norm) of the vector</returns>
|
|
|
|
private static float Cnrm2Vector(Complex32[] a, int rowStart) |
|
|
|
static float Cnrm2Vector(Complex32[] a, int rowStart) |
|
|
|
{ |
|
|
|
var s = 0.0f; |
|
|
|
for (var i = rowStart; i < a.Length; i++) |
|
|
|
{ |
|
|
|
s += a[i].Magnitude * a[i].Magnitude; |
|
|
|
s += a[i].Magnitude*a[i].Magnitude; |
|
|
|
} |
|
|
|
|
|
|
|
return (float)Math.Sqrt(s); |
|
|
|
return (float) Math.Sqrt(s); |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
@ -763,12 +764,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="columnB">Index of column B</param>
|
|
|
|
/// <param name="rowStart">Starting row index</param>
|
|
|
|
/// <returns>Dot product value</returns>
|
|
|
|
private static Complex32 Cdotc(Matrix<Complex32> a, int rowCount, int columnA, int columnB, int rowStart) |
|
|
|
static Complex32 Cdotc(Matrix<Complex32> a, int rowCount, int columnA, int columnB, int rowStart) |
|
|
|
{ |
|
|
|
var z = Complex32.Zero; |
|
|
|
for (var i = rowStart; i < rowCount; i++) |
|
|
|
{ |
|
|
|
z += a.At(i, columnA).Conjugate() * a.At(i, columnB); |
|
|
|
z += a.At(i, columnA).Conjugate()*a.At(i, columnB); |
|
|
|
} |
|
|
|
|
|
|
|
return z; |
|
|
|
@ -784,17 +785,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="columnB">Index of column B</param>
|
|
|
|
/// <param name="c">scalar cos value</param>
|
|
|
|
/// <param name="s">scalar sin value</param>
|
|
|
|
private static void Csrot(Matrix<Complex32> a, int rowCount, int columnA, int columnB, float c, float s) |
|
|
|
static void Csrot(Matrix<Complex32> a, int rowCount, int columnA, int columnB, float c, float s) |
|
|
|
{ |
|
|
|
for (var i = 0; i < rowCount; i++) |
|
|
|
{ |
|
|
|
var z = (c * a.At(i, columnA)) + (s * a.At(i, columnB)); |
|
|
|
var tmp = (c * a.At(i, columnB)) - (s * a.At(i, columnA)); |
|
|
|
var z = (c*a.At(i, columnA)) + (s*a.At(i, columnB)); |
|
|
|
var tmp = (c*a.At(i, columnB)) - (s*a.At(i, columnA)); |
|
|
|
a.At(i, columnB, tmp); |
|
|
|
a.At(i, columnA, z); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Solves a system of linear equations, <b>AX = B</b>, with A SVD factorized.
|
|
|
|
/// </summary>
|
|
|
|
@ -802,18 +803,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>X</b>.</param>
|
|
|
|
public override void Solve(Matrix<Complex32> input, Matrix<Complex32> result) |
|
|
|
{ |
|
|
|
// Check for proper arguments.
|
|
|
|
if (input == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("input"); |
|
|
|
} |
|
|
|
|
|
|
|
if (result == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("result"); |
|
|
|
} |
|
|
|
|
|
|
|
if (!ComputeVectors) |
|
|
|
if (!VectorsComputed) |
|
|
|
{ |
|
|
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
|
|
|
} |
|
|
|
@ -850,7 +840,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
{ |
|
|
|
for (var i = 0; i < U.RowCount; i++) |
|
|
|
{ |
|
|
|
value += U.At(i, j).Conjugate() * input.At(i, k); |
|
|
|
value += U.At(i, j).Conjugate()*input.At(i, k); |
|
|
|
} |
|
|
|
|
|
|
|
value /= S[j]; |
|
|
|
@ -864,7 +854,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
var value = Complex32.Zero; |
|
|
|
for (var i = 0; i < VT.ColumnCount; i++) |
|
|
|
{ |
|
|
|
value += VT.At(i, j).Conjugate() * tmp[i]; |
|
|
|
value += VT.At(i, j).Conjugate()*tmp[i]; |
|
|
|
} |
|
|
|
|
|
|
|
result.At(j, k, value); |
|
|
|
@ -879,17 +869,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
/// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
|
|
|
|
public override void Solve(Vector<Complex32> input, Vector<Complex32> result) |
|
|
|
{ |
|
|
|
if (input == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("input"); |
|
|
|
} |
|
|
|
|
|
|
|
if (result == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("result"); |
|
|
|
} |
|
|
|
|
|
|
|
if (!ComputeVectors) |
|
|
|
if (!VectorsComputed) |
|
|
|
{ |
|
|
|
throw new InvalidOperationException(Resources.SingularVectorsNotComputed); |
|
|
|
} |
|
|
|
@ -916,7 +896,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
{ |
|
|
|
for (var i = 0; i < U.RowCount; i++) |
|
|
|
{ |
|
|
|
value += U.At(i, j).Conjugate() * input[i]; |
|
|
|
value += U.At(i, j).Conjugate()*input[i]; |
|
|
|
} |
|
|
|
|
|
|
|
value /= S[j]; |
|
|
|
@ -930,7 +910,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization |
|
|
|
var value = Complex32.Zero; |
|
|
|
for (var i = 0; i < VT.ColumnCount; i++) |
|
|
|
{ |
|
|
|
value += VT.At(i, j).Conjugate() * tmp[i]; |
|
|
|
value += VT.At(i, j).Conjugate()*tmp[i]; |
|
|
|
} |
|
|
|
|
|
|
|
result[j] = value; |
|
|
|
|
Christoph Ruegg
13 years ago