mirror of
https://github.com/OpenFOAM/ThirdParty-6.git
synced 2025-12-08 06:57:43 +00:00
ParaView-5.0.1: Added the source-tree to ThirdParty-dev and patched as described in the README file
Resolves bug-report http://bugs.openfoam.org/view.php?id=2098
This commit is contained in:
Henry Weller
611
ParaView-5.0.1/VTK/ThirdParty/alglib/svd.cpp
vendored
Normal file
611
ParaView-5.0.1/VTK/ThirdParty/alglib/svd.cpp
vendored
Normal file
@ -0,0 +1,611 @@
|
||||
/*************************************************************************
|
||||
Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
|
||||
|
||||
Redistribution and use in source and binary forms, with or without
|
||||
modification, are permitted provided that the following conditions are
|
||||
met:
|
||||
|
||||
- Redistributions of source code must retain the above copyright
|
||||
notice, this list of conditions and the following disclaimer.
|
||||
|
||||
- Redistributions in binary form must reproduce the above copyright
|
||||
notice, this list of conditions and the following disclaimer listed
|
||||
in this license in the documentation and/or other materials
|
||||
provided with the distribution.
|
||||
|
||||
- Neither the name of the copyright holders nor the names of its
|
||||
contributors may be used to endorse or promote products derived from
|
||||
this software without specific prior written permission.
|
||||
|
||||
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
||||
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
||||
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
||||
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
||||
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
||||
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
||||
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
||||
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
||||
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
||||
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
||||
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||||
*************************************************************************/
|
||||
|
||||
#include "alglib/svd.h"
|
||||
|
||||
/*************************************************************************
|
||||
Singular value decomposition of a rectangular matrix.
|
||||
|
||||
The algorithm calculates the singular value decomposition of a matrix of
|
||||
size MxN: A = U * S * V^T
|
||||
|
||||
The algorithm finds the singular values and, optionally, matrices U and V^T.
|
||||
The algorithm can find both first min(M,N) columns of matrix U and rows of
|
||||
matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
|
||||
and NxN respectively).
|
||||
|
||||
Take into account that the subroutine does not return matrix V but V^T.
|
||||
|
||||
Input parameters:
|
||||
A - matrix to be decomposed.
|
||||
Array whose indexes range within [0..M-1, 0..N-1].
|
||||
M - number of rows in matrix A.
|
||||
N - number of columns in matrix A.
|
||||
UNeeded - 0, 1 or 2. See the description of the parameter U.
|
||||
VTNeeded - 0, 1 or 2. See the description of the parameter VT.
|
||||
AdditionalMemory -
|
||||
If the parameter:
|
||||
* equals 0, the algorithm doesn<73>t use additional
|
||||
memory (lower requirements, lower performance).
|
||||
* equals 1, the algorithm uses additional
|
||||
memory of size min(M,N)*min(M,N) of real numbers.
|
||||
It often speeds up the algorithm.
|
||||
* equals 2, the algorithm uses additional
|
||||
memory of size M*min(M,N) of real numbers.
|
||||
It allows to get a maximum performance.
|
||||
The recommended value of the parameter is 2.
|
||||
|
||||
Output parameters:
|
||||
W - contains singular values in descending order.
|
||||
U - if UNeeded=0, U isn't changed, the left singular vectors
|
||||
are not calculated.
|
||||
if Uneeded=1, U contains left singular vectors (first
|
||||
min(M,N) columns of matrix U). Array whose indexes range
|
||||
within [0..M-1, 0..Min(M,N)-1].
|
||||
if UNeeded=2, U contains matrix U wholly. Array whose
|
||||
indexes range within [0..M-1, 0..M-1].
|
||||
VT - if VTNeeded=0, VT isn<73>t changed, the right singular vectors
|
||||
are not calculated.
|
||||
if VTNeeded=1, VT contains right singular vectors (first
|
||||
min(M,N) rows of matrix V^T). Array whose indexes range
|
||||
within [0..min(M,N)-1, 0..N-1].
|
||||
if VTNeeded=2, VT contains matrix V^T wholly. Array whose
|
||||
indexes range within [0..N-1, 0..N-1].
|
||||
|
||||
-- ALGLIB --
|
||||
Copyright 2005 by Bochkanov Sergey
|
||||
*************************************************************************/
|
||||
bool rmatrixsvd(ap::real_2d_array a,
|
||||
int m,
|
||||
int n,
|
||||
int uneeded,
|
||||
int vtneeded,
|
||||
int additionalmemory,
|
||||
ap::real_1d_array& w,
|
||||
ap::real_2d_array& u,
|
||||
ap::real_2d_array& vt)
|
||||
{
|
||||
bool result;
|
||||
ap::real_1d_array tauq;
|
||||
ap::real_1d_array taup;
|
||||
ap::real_1d_array tau;
|
||||
ap::real_1d_array e;
|
||||
ap::real_1d_array work;
|
||||
ap::real_2d_array t2;
|
||||
bool isupper;
|
||||
int minmn;
|
||||
int ncu;
|
||||
int nrvt;
|
||||
int nru;
|
||||
int ncvt;
|
||||
int i;
|
||||
int j;
|
||||
|
||||
result = true;
|
||||
if( m==0||n==0 )
|
||||
{
|
||||
return result;
|
||||
}
|
||||
ap::ap_error::make_assertion(uneeded>=0&&uneeded<=2, "SVDDecomposition: wrong parameters!");
|
||||
ap::ap_error::make_assertion(vtneeded>=0&&vtneeded<=2, "SVDDecomposition: wrong parameters!");
|
||||
ap::ap_error::make_assertion(additionalmemory>=0&&additionalmemory<=2, "SVDDecomposition: wrong parameters!");
|
||||
|
||||
//
|
||||
// initialize
|
||||
//
|
||||
minmn = ap::minint(m, n);
|
||||
w.setbounds(1, minmn);
|
||||
ncu = 0;
|
||||
nru = 0;
|
||||
if( uneeded==1 )
|
||||
{
|
||||
nru = m;
|
||||
ncu = minmn;
|
||||
u.setbounds(0, nru-1, 0, ncu-1);
|
||||
}
|
||||
if( uneeded==2 )
|
||||
{
|
||||
nru = m;
|
||||
ncu = m;
|
||||
u.setbounds(0, nru-1, 0, ncu-1);
|
||||
}
|
||||
nrvt = 0;
|
||||
ncvt = 0;
|
||||
if( vtneeded==1 )
|
||||
{
|
||||
nrvt = minmn;
|
||||
ncvt = n;
|
||||
vt.setbounds(0, nrvt-1, 0, ncvt-1);
|
||||
}
|
||||
if( vtneeded==2 )
|
||||
{
|
||||
nrvt = n;
|
||||
ncvt = n;
|
||||
vt.setbounds(0, nrvt-1, 0, ncvt-1);
|
||||
}
|
||||
|
||||
//
|
||||
// M much larger than N
|
||||
// Use bidiagonal reduction with QR-decomposition
|
||||
//
|
||||
if( m>1.6*n )
|
||||
{
|
||||
if( uneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// No left singular vectors to be computed
|
||||
//
|
||||
rmatrixqr(a, m, n, tau);
|
||||
for(i = 0; i <= n-1; i++)
|
||||
{
|
||||
for(j = 0; j <= i-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
rmatrixbd(a, n, n, tauq, taup);
|
||||
rmatrixbdunpackpt(a, n, n, taup, nrvt, vt);
|
||||
rmatrixbdunpackdiagonals(a, n, n, isupper, w, e);
|
||||
result = rmatrixbdsvd(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
|
||||
return result;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Left singular vectors (may be full matrix U) to be computed
|
||||
//
|
||||
rmatrixqr(a, m, n, tau);
|
||||
rmatrixqrunpackq(a, m, n, tau, ncu, u);
|
||||
for(i = 0; i <= n-1; i++)
|
||||
{
|
||||
for(j = 0; j <= i-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
rmatrixbd(a, n, n, tauq, taup);
|
||||
rmatrixbdunpackpt(a, n, n, taup, nrvt, vt);
|
||||
rmatrixbdunpackdiagonals(a, n, n, isupper, w, e);
|
||||
if( additionalmemory<1 )
|
||||
{
|
||||
|
||||
//
|
||||
// No additional memory can be used
|
||||
//
|
||||
rmatrixbdmultiplybyq(a, n, n, tauq, u, m, n, true, false);
|
||||
result = rmatrixbdsvd(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Large U. Transforming intermediate matrix T2
|
||||
//
|
||||
work.setbounds(1, ap::maxint(m, n));
|
||||
rmatrixbdunpackq(a, n, n, tauq, n, t2);
|
||||
copymatrix(u, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
|
||||
inplacetranspose(t2, 0, n-1, 0, n-1, work);
|
||||
result = rmatrixbdsvd(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
|
||||
matrixmatrixmultiply(a, 0, m-1, 0, n-1, false, t2, 0, n-1, 0, n-1, true, 1.0, u, 0, m-1, 0, n-1, 0.0, work);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
// N much larger than M
|
||||
// Use bidiagonal reduction with LQ-decomposition
|
||||
//
|
||||
if( n>1.6*m )
|
||||
{
|
||||
if( vtneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// No right singular vectors to be computed
|
||||
//
|
||||
rmatrixlq(a, m, n, tau);
|
||||
for(i = 0; i <= m-1; i++)
|
||||
{
|
||||
for(j = i+1; j <= m-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
rmatrixbd(a, m, m, tauq, taup);
|
||||
rmatrixbdunpackq(a, m, m, tauq, ncu, u);
|
||||
rmatrixbdunpackdiagonals(a, m, m, isupper, w, e);
|
||||
work.setbounds(1, m);
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
result = rmatrixbdsvd(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
return result;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Right singular vectors (may be full matrix VT) to be computed
|
||||
//
|
||||
rmatrixlq(a, m, n, tau);
|
||||
rmatrixlqunpackq(a, m, n, tau, nrvt, vt);
|
||||
for(i = 0; i <= m-1; i++)
|
||||
{
|
||||
for(j = i+1; j <= m-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
rmatrixbd(a, m, m, tauq, taup);
|
||||
rmatrixbdunpackq(a, m, m, tauq, ncu, u);
|
||||
rmatrixbdunpackdiagonals(a, m, m, isupper, w, e);
|
||||
work.setbounds(1, ap::maxint(m, n));
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
if( additionalmemory<1 )
|
||||
{
|
||||
|
||||
//
|
||||
// No additional memory available
|
||||
//
|
||||
rmatrixbdmultiplybyp(a, m, m, taup, vt, m, n, false, true);
|
||||
result = rmatrixbdsvd(w, e, m, isupper, false, a, 0, u, nru, vt, n);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Large VT. Transforming intermediate matrix T2
|
||||
//
|
||||
rmatrixbdunpackpt(a, m, m, taup, m, t2);
|
||||
result = rmatrixbdsvd(w, e, m, isupper, false, a, 0, u, nru, t2, m);
|
||||
copymatrix(vt, 0, m-1, 0, n-1, a, 0, m-1, 0, n-1);
|
||||
matrixmatrixmultiply(t2, 0, m-1, 0, m-1, false, a, 0, m-1, 0, n-1, false, 1.0, vt, 0, m-1, 0, n-1, 0.0, work);
|
||||
}
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
// M<=N
|
||||
// We can use inplace transposition of U to get rid of columnwise operations
|
||||
//
|
||||
if( m<=n )
|
||||
{
|
||||
rmatrixbd(a, m, n, tauq, taup);
|
||||
rmatrixbdunpackq(a, m, n, tauq, ncu, u);
|
||||
rmatrixbdunpackpt(a, m, n, taup, nrvt, vt);
|
||||
rmatrixbdunpackdiagonals(a, m, n, isupper, w, e);
|
||||
work.setbounds(1, m);
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
result = rmatrixbdsvd(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
|
||||
inplacetranspose(u, 0, nru-1, 0, ncu-1, work);
|
||||
return result;
|
||||
}
|
||||
|
||||
//
|
||||
// Simple bidiagonal reduction
|
||||
//
|
||||
rmatrixbd(a, m, n, tauq, taup);
|
||||
rmatrixbdunpackq(a, m, n, tauq, ncu, u);
|
||||
rmatrixbdunpackpt(a, m, n, taup, nrvt, vt);
|
||||
rmatrixbdunpackdiagonals(a, m, n, isupper, w, e);
|
||||
if( additionalmemory<2||uneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// We cant use additional memory or there is no need in such operations
|
||||
//
|
||||
result = rmatrixbdsvd(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// We can use additional memory
|
||||
//
|
||||
t2.setbounds(0, minmn-1, 0, m-1);
|
||||
copyandtranspose(u, 0, m-1, 0, minmn-1, t2, 0, minmn-1, 0, m-1);
|
||||
result = rmatrixbdsvd(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
|
||||
copyandtranspose(t2, 0, minmn-1, 0, m-1, u, 0, m-1, 0, minmn-1);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
/*************************************************************************
|
||||
Obsolete 1-based subroutine.
|
||||
See RMatrixSVD for 0-based replacement.
|
||||
*************************************************************************/
|
||||
bool svddecomposition(ap::real_2d_array a,
|
||||
int m,
|
||||
int n,
|
||||
int uneeded,
|
||||
int vtneeded,
|
||||
int additionalmemory,
|
||||
ap::real_1d_array& w,
|
||||
ap::real_2d_array& u,
|
||||
ap::real_2d_array& vt)
|
||||
{
|
||||
bool result;
|
||||
ap::real_1d_array tauq;
|
||||
ap::real_1d_array taup;
|
||||
ap::real_1d_array tau;
|
||||
ap::real_1d_array e;
|
||||
ap::real_1d_array work;
|
||||
ap::real_2d_array t2;
|
||||
bool isupper;
|
||||
int minmn;
|
||||
int ncu;
|
||||
int nrvt;
|
||||
int nru;
|
||||
int ncvt;
|
||||
int i;
|
||||
int j;
|
||||
|
||||
result = true;
|
||||
if( m==0||n==0 )
|
||||
{
|
||||
return result;
|
||||
}
|
||||
ap::ap_error::make_assertion(uneeded>=0&&uneeded<=2, "SVDDecomposition: wrong parameters!");
|
||||
ap::ap_error::make_assertion(vtneeded>=0&&vtneeded<=2, "SVDDecomposition: wrong parameters!");
|
||||
ap::ap_error::make_assertion(additionalmemory>=0&&additionalmemory<=2, "SVDDecomposition: wrong parameters!");
|
||||
|
||||
//
|
||||
// initialize
|
||||
//
|
||||
minmn = ap::minint(m, n);
|
||||
w.setbounds(1, minmn);
|
||||
ncu = 0;
|
||||
nru = 0;
|
||||
if( uneeded==1 )
|
||||
{
|
||||
nru = m;
|
||||
ncu = minmn;
|
||||
u.setbounds(1, nru, 1, ncu);
|
||||
}
|
||||
if( uneeded==2 )
|
||||
{
|
||||
nru = m;
|
||||
ncu = m;
|
||||
u.setbounds(1, nru, 1, ncu);
|
||||
}
|
||||
nrvt = 0;
|
||||
ncvt = 0;
|
||||
if( vtneeded==1 )
|
||||
{
|
||||
nrvt = minmn;
|
||||
ncvt = n;
|
||||
vt.setbounds(1, nrvt, 1, ncvt);
|
||||
}
|
||||
if( vtneeded==2 )
|
||||
{
|
||||
nrvt = n;
|
||||
ncvt = n;
|
||||
vt.setbounds(1, nrvt, 1, ncvt);
|
||||
}
|
||||
|
||||
//
|
||||
// M much larger than N
|
||||
// Use bidiagonal reduction with QR-decomposition
|
||||
//
|
||||
if( m>1.6*n )
|
||||
{
|
||||
if( uneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// No left singular vectors to be computed
|
||||
//
|
||||
qrdecomposition(a, m, n, tau);
|
||||
for(i = 2; i <= n; i++)
|
||||
{
|
||||
for(j = 1; j <= i-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
tobidiagonal(a, n, n, tauq, taup);
|
||||
unpackptfrombidiagonal(a, n, n, taup, nrvt, vt);
|
||||
unpackdiagonalsfrombidiagonal(a, n, n, isupper, w, e);
|
||||
result = bidiagonalsvddecomposition(w, e, n, isupper, false, u, 0, a, 0, vt, ncvt);
|
||||
return result;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Left singular vectors (may be full matrix U) to be computed
|
||||
//
|
||||
qrdecomposition(a, m, n, tau);
|
||||
unpackqfromqr(a, m, n, tau, ncu, u);
|
||||
for(i = 2; i <= n; i++)
|
||||
{
|
||||
for(j = 1; j <= i-1; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
tobidiagonal(a, n, n, tauq, taup);
|
||||
unpackptfrombidiagonal(a, n, n, taup, nrvt, vt);
|
||||
unpackdiagonalsfrombidiagonal(a, n, n, isupper, w, e);
|
||||
if( additionalmemory<1 )
|
||||
{
|
||||
|
||||
//
|
||||
// No additional memory can be used
|
||||
//
|
||||
multiplybyqfrombidiagonal(a, n, n, tauq, u, m, n, true, false);
|
||||
result = bidiagonalsvddecomposition(w, e, n, isupper, false, u, m, a, 0, vt, ncvt);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Large U. Transforming intermediate matrix T2
|
||||
//
|
||||
work.setbounds(1, ap::maxint(m, n));
|
||||
unpackqfrombidiagonal(a, n, n, tauq, n, t2);
|
||||
copymatrix(u, 1, m, 1, n, a, 1, m, 1, n);
|
||||
inplacetranspose(t2, 1, n, 1, n, work);
|
||||
result = bidiagonalsvddecomposition(w, e, n, isupper, false, u, 0, t2, n, vt, ncvt);
|
||||
matrixmatrixmultiply(a, 1, m, 1, n, false, t2, 1, n, 1, n, true, 1.0, u, 1, m, 1, n, 0.0, work);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
// N much larger than M
|
||||
// Use bidiagonal reduction with LQ-decomposition
|
||||
//
|
||||
if( n>1.6*m )
|
||||
{
|
||||
if( vtneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// No right singular vectors to be computed
|
||||
//
|
||||
lqdecomposition(a, m, n, tau);
|
||||
for(i = 1; i <= m-1; i++)
|
||||
{
|
||||
for(j = i+1; j <= m; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
tobidiagonal(a, m, m, tauq, taup);
|
||||
unpackqfrombidiagonal(a, m, m, tauq, ncu, u);
|
||||
unpackdiagonalsfrombidiagonal(a, m, m, isupper, w, e);
|
||||
work.setbounds(1, m);
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
result = bidiagonalsvddecomposition(w, e, m, isupper, false, a, 0, u, nru, vt, 0);
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
return result;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Right singular vectors (may be full matrix VT) to be computed
|
||||
//
|
||||
lqdecomposition(a, m, n, tau);
|
||||
unpackqfromlq(a, m, n, tau, nrvt, vt);
|
||||
for(i = 1; i <= m-1; i++)
|
||||
{
|
||||
for(j = i+1; j <= m; j++)
|
||||
{
|
||||
a(i,j) = 0;
|
||||
}
|
||||
}
|
||||
tobidiagonal(a, m, m, tauq, taup);
|
||||
unpackqfrombidiagonal(a, m, m, tauq, ncu, u);
|
||||
unpackdiagonalsfrombidiagonal(a, m, m, isupper, w, e);
|
||||
work.setbounds(1, ap::maxint(m, n));
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
if( additionalmemory<1 )
|
||||
{
|
||||
|
||||
//
|
||||
// No additional memory available
|
||||
//
|
||||
multiplybypfrombidiagonal(a, m, m, taup, vt, m, n, false, true);
|
||||
result = bidiagonalsvddecomposition(w, e, m, isupper, false, a, 0, u, nru, vt, n);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// Large VT. Transforming intermediate matrix T2
|
||||
//
|
||||
unpackptfrombidiagonal(a, m, m, taup, m, t2);
|
||||
result = bidiagonalsvddecomposition(w, e, m, isupper, false, a, 0, u, nru, t2, m);
|
||||
copymatrix(vt, 1, m, 1, n, a, 1, m, 1, n);
|
||||
matrixmatrixmultiply(t2, 1, m, 1, m, false, a, 1, m, 1, n, false, 1.0, vt, 1, m, 1, n, 0.0, work);
|
||||
}
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
//
|
||||
// M<=N
|
||||
// We can use inplace transposition of U to get rid of columnwise operations
|
||||
//
|
||||
if( m<=n )
|
||||
{
|
||||
tobidiagonal(a, m, n, tauq, taup);
|
||||
unpackqfrombidiagonal(a, m, n, tauq, ncu, u);
|
||||
unpackptfrombidiagonal(a, m, n, taup, nrvt, vt);
|
||||
unpackdiagonalsfrombidiagonal(a, m, n, isupper, w, e);
|
||||
work.setbounds(1, m);
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
result = bidiagonalsvddecomposition(w, e, minmn, isupper, false, a, 0, u, nru, vt, ncvt);
|
||||
inplacetranspose(u, 1, nru, 1, ncu, work);
|
||||
return result;
|
||||
}
|
||||
|
||||
//
|
||||
// Simple bidiagonal reduction
|
||||
//
|
||||
tobidiagonal(a, m, n, tauq, taup);
|
||||
unpackqfrombidiagonal(a, m, n, tauq, ncu, u);
|
||||
unpackptfrombidiagonal(a, m, n, taup, nrvt, vt);
|
||||
unpackdiagonalsfrombidiagonal(a, m, n, isupper, w, e);
|
||||
if( additionalmemory<2||uneeded==0 )
|
||||
{
|
||||
|
||||
//
|
||||
// We cant use additional memory or there is no need in such operations
|
||||
//
|
||||
result = bidiagonalsvddecomposition(w, e, minmn, isupper, false, u, nru, a, 0, vt, ncvt);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
||||
//
|
||||
// We can use additional memory
|
||||
//
|
||||
t2.setbounds(1, minmn, 1, m);
|
||||
copyandtranspose(u, 1, m, 1, minmn, t2, 1, minmn, 1, m);
|
||||
result = bidiagonalsvddecomposition(w, e, minmn, isupper, false, u, 0, t2, m, vt, ncvt);
|
||||
copyandtranspose(t2, 1, minmn, 1, m, u, 1, m, 1, minmn);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user