Context Navigation

source: trunk/sources/ALGLIB/svd.cs @ 2563

Visit:

Last change on this file since 2563 was 2563, checked in by gkronber, 15 years ago
Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 26.2 KB

Rev	Line
[2154]	1	/*************************************************************************
	2	Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
	3
[2430]	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
[2154]	9
[2430]	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
[2154]	14
[2430]	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
[2154]	17
[2430]	18	>>> END OF LICENSE >>>
[2154]	19	*************************************************************************/
	20
	21	using System;
	22
[2430]	23	namespace alglib
[2154]	24	{
[2430]	25	public class svd
	26	{
	27	/*************************************************************************
	28	Singular value decomposition of a rectangular matrix.
[2154]	29
[2430]	30	The algorithm calculates the singular value decomposition of a matrix of
	31	size MxN: A = U * S * V^T
[2154]	32
[2430]	33	The algorithm finds the singular values and, optionally, matrices U and V^T.
	34	The algorithm can find both first min(M,N) columns of matrix U and rows of
	35	matrix V^T (singular vectors), and matrices U and V^T wholly (of sizes MxM
	36	and NxN respectively).
[2154]	37
[2430]	38	Take into account that the subroutine does not return matrix V but V^T.
[2154]	39
[2430]	40	Input parameters:
	41	A - matrix to be decomposed.
	42	Array whose indexes range within [0..M-1, 0..N-1].
	43	M - number of rows in matrix A.
	44	N - number of columns in matrix A.
	45	UNeeded - 0, 1 or 2. See the description of the parameter U.
	46	VTNeeded - 0, 1 or 2. See the description of the parameter VT.
	47	AdditionalMemory -
	48	If the parameter:
	49	* equals 0, the algorithm doesnt use additional
	50	memory (lower requirements, lower performance).
	51	* equals 1, the algorithm uses additional
	52	memory of size min(M,N)*min(M,N) of real numbers.
	53	It often speeds up the algorithm.
	54	* equals 2, the algorithm uses additional
	55	memory of size M*min(M,N) of real numbers.
	56	It allows to get a maximum performance.
	57	The recommended value of the parameter is 2.
[2154]	58
[2430]	59	Output parameters:
	60	W - contains singular values in descending order.
	61	U - if UNeeded=0, U isn't changed, the left singular vectors
	62	are not calculated.
	63	if Uneeded=1, U contains left singular vectors (first
	64	min(M,N) columns of matrix U). Array whose indexes range
	65	within [0..M-1, 0..Min(M,N)-1].
	66	if UNeeded=2, U contains matrix U wholly. Array whose
	67	indexes range within [0..M-1, 0..M-1].
	68	VT - if VTNeeded=0, VT isnt changed, the right singular vectors
	69	are not calculated.
	70	if VTNeeded=1, VT contains right singular vectors (first
	71	min(M,N) rows of matrix V^T). Array whose indexes range
	72	within [0..min(M,N)-1, 0..N-1].
	73	if VTNeeded=2, VT contains matrix V^T wholly. Array whose
	74	indexes range within [0..N-1, 0..N-1].
[2154]	75
[2430]	76	-- ALGLIB --
	77	Copyright 2005 by Bochkanov Sergey
	78	*************************************************************************/
	79	public static bool rmatrixsvd(double[,] a,
	80	int m,
	81	int n,
	82	int uneeded,
	83	int vtneeded,
	84	int additionalmemory,
	85	ref double[] w,
	86	ref double[,] u,
	87	ref double[,] vt)
	88	{
	89	bool result = new bool();
	90	double[] tauq = new double[0];
	91	double[] taup = new double[0];
	92	double[] tau = new double[0];
	93	double[] e = new double[0];
	94	double[] work = new double[0];
	95	double[,] t2 = new double[0,0];
	96	bool isupper = new bool();
	97	int minmn = 0;
	98	int ncu = 0;
	99	int nrvt = 0;
	100	int nru = 0;
	101	int ncvt = 0;
	102	int i = 0;
	103	int j = 0;
	104	int im1 = 0;
	105	double sm = 0;
[2154]	106
[2430]	107	a = (double[,])a.Clone();
[2154]	108
[2430]	109	result = true;
	110	if( m==0 \| n==0 )
[2154]	111	{
	112	return result;
	113	}
[2430]	114	System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
	115	System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
	116	System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
	117
	118	//
	119	// initialize
	120	//
	121	minmn = Math.Min(m, n);
	122	w = new double[minmn+1];
	123	ncu = 0;
	124	nru = 0;
	125	if( uneeded==1 )
[2154]	126	{
[2430]	127	nru = m;
	128	ncu = minmn;
	129	u = new double[nru-1+1, ncu-1+1];
	130	}
	131	if( uneeded==2 )
	132	{
	133	nru = m;
	134	ncu = m;
	135	u = new double[nru-1+1, ncu-1+1];
	136	}
	137	nrvt = 0;
	138	ncvt = 0;
	139	if( vtneeded==1 )
	140	{
	141	nrvt = minmn;
	142	ncvt = n;
	143	vt = new double[nrvt-1+1, ncvt-1+1];
	144	}
	145	if( vtneeded==2 )
	146	{
	147	nrvt = n;
	148	ncvt = n;
	149	vt = new double[nrvt-1+1, ncvt-1+1];
	150	}
	151
	152	//
	153	// M much larger than N
	154	// Use bidiagonal reduction with QR-decomposition
	155	//
[2563]	156	if( (double)(m)>(double)(1.6*n) )
[2430]	157	{
	158	if( uneeded==0 )
[2154]	159	{
	160
	161	//
[2430]	162	// No left singular vectors to be computed
[2154]	163	//
[2430]	164	qr.rmatrixqr(ref a, m, n, ref tau);
	165	for(i=0; i<=n-1; i++)
	166	{
	167	for(j=0; j<=i-1; j++)
	168	{
	169	a[i,j] = 0;
	170	}
	171	}
	172	bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
	173	bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
	174	bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
	175	result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
	176	return result;
[2154]	177	}
	178	else
	179	{
	180
	181	//
[2430]	182	// Left singular vectors (may be full matrix U) to be computed
[2154]	183	//
[2430]	184	qr.rmatrixqr(ref a, m, n, ref tau);
	185	qr.rmatrixqrunpackq(ref a, m, n, ref tau, ncu, ref u);
	186	for(i=0; i<=n-1; i++)
[2154]	187	{
[2430]	188	for(j=0; j<=i-1; j++)
	189	{
	190	a[i,j] = 0;
	191	}
[2154]	192	}
[2430]	193	bidiagonal.rmatrixbd(ref a, n, n, ref tauq, ref taup);
	194	bidiagonal.rmatrixbdunpackpt(ref a, n, n, ref taup, nrvt, ref vt);
	195	bidiagonal.rmatrixbdunpackdiagonals(ref a, n, n, ref isupper, ref w, ref e);
	196	if( additionalmemory<1 )
	197	{
	198
	199	//
	200	// No additional memory can be used
	201	//
	202	bidiagonal.rmatrixbdmultiplybyq(ref a, n, n, ref tauq, ref u, m, n, true, false);
	203	result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
	204	}
	205	else
	206	{
	207
	208	//
	209	// Large U. Transforming intermediate matrix T2
	210	//
	211	work = new double[Math.Max(m, n)+1];
	212	bidiagonal.rmatrixbdunpackq(ref a, n, n, ref tauq, n, ref t2);
	213	blas.copymatrix(ref u, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
	214	blas.inplacetranspose(ref t2, 0, n-1, 0, n-1, ref work);
	215	result = bdsvd.rmatrixbdsvd(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
	216	blas.matrixmatrixmultiply(ref a, 0, m-1, 0, n-1, false, ref t2, 0, n-1, 0, n-1, true, 1.0, ref u, 0, m-1, 0, n-1, 0.0, ref work);
	217	}
	218	return result;
[2154]	219	}
	220	}
[2430]	221
	222	//
	223	// N much larger than M
	224	// Use bidiagonal reduction with LQ-decomposition
	225	//
[2563]	226	if( (double)(n)>(double)(1.6*m) )
[2154]	227	{
[2430]	228	if( vtneeded==0 )
[2154]	229	{
	230
	231	//
[2430]	232	// No right singular vectors to be computed
[2154]	233	//
[2430]	234	lq.rmatrixlq(ref a, m, n, ref tau);
	235	for(i=0; i<=m-1; i++)
	236	{
	237	for(j=i+1; j<=m-1; j++)
	238	{
	239	a[i,j] = 0;
	240	}
	241	}
	242	bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
	243	bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
	244	bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
	245	work = new double[m+1];
	246	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
	247	result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
	248	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
	249	return result;
[2154]	250	}
	251	else
	252	{
	253
	254	//
[2430]	255	// Right singular vectors (may be full matrix VT) to be computed
[2154]	256	//
[2430]	257	lq.rmatrixlq(ref a, m, n, ref tau);
	258	lq.rmatrixlqunpackq(ref a, m, n, ref tau, nrvt, ref vt);
	259	for(i=0; i<=m-1; i++)
	260	{
	261	for(j=i+1; j<=m-1; j++)
	262	{
	263	a[i,j] = 0;
	264	}
	265	}
	266	bidiagonal.rmatrixbd(ref a, m, m, ref tauq, ref taup);
	267	bidiagonal.rmatrixbdunpackq(ref a, m, m, ref tauq, ncu, ref u);
	268	bidiagonal.rmatrixbdunpackdiagonals(ref a, m, m, ref isupper, ref w, ref e);
	269	work = new double[Math.Max(m, n)+1];
	270	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
	271	if( additionalmemory<1 )
	272	{
	273
	274	//
	275	// No additional memory available
	276	//
	277	bidiagonal.rmatrixbdmultiplybyp(ref a, m, m, ref taup, ref vt, m, n, false, true);
	278	result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
	279	}
	280	else
	281	{
	282
	283	//
	284	// Large VT. Transforming intermediate matrix T2
	285	//
	286	bidiagonal.rmatrixbdunpackpt(ref a, m, m, ref taup, m, ref t2);
	287	result = bdsvd.rmatrixbdsvd(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
	288	blas.copymatrix(ref vt, 0, m-1, 0, n-1, ref a, 0, m-1, 0, n-1);
	289	blas.matrixmatrixmultiply(ref t2, 0, m-1, 0, m-1, false, ref a, 0, m-1, 0, n-1, false, 1.0, ref vt, 0, m-1, 0, n-1, 0.0, ref work);
	290	}
	291	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
	292	return result;
[2154]	293	}
[2430]	294	}
	295
	296	//
	297	// M<=N
	298	// We can use inplace transposition of U to get rid of columnwise operations
	299	//
	300	if( m<=n )
	301	{
	302	bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
	303	bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
	304	bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
	305	bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
	306	work = new double[m+1];
[2154]	307	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
[2430]	308	result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
	309	blas.inplacetranspose(ref u, 0, nru-1, 0, ncu-1, ref work);
[2154]	310	return result;
	311	}
[2430]	312
	313	//
	314	// Simple bidiagonal reduction
	315	//
[2154]	316	bidiagonal.rmatrixbd(ref a, m, n, ref tauq, ref taup);
	317	bidiagonal.rmatrixbdunpackq(ref a, m, n, ref tauq, ncu, ref u);
	318	bidiagonal.rmatrixbdunpackpt(ref a, m, n, ref taup, nrvt, ref vt);
	319	bidiagonal.rmatrixbdunpackdiagonals(ref a, m, n, ref isupper, ref w, ref e);
[2430]	320	if( additionalmemory<2 \| uneeded==0 )
	321	{
	322
	323	//
	324	// We cant use additional memory or there is no need in such operations
	325	//
	326	result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
	327	}
	328	else
	329	{
	330
	331	//
	332	// We can use additional memory
	333	//
	334	t2 = new double[minmn-1+1, m-1+1];
	335	blas.copyandtranspose(ref u, 0, m-1, 0, minmn-1, ref t2, 0, minmn-1, 0, m-1);
	336	result = bdsvd.rmatrixbdsvd(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
	337	blas.copyandtranspose(ref t2, 0, minmn-1, 0, m-1, ref u, 0, m-1, 0, minmn-1);
	338	}
[2154]	339	return result;
	340	}
[2430]	341
	342
	343	public static bool svddecomposition(double[,] a,
	344	int m,
	345	int n,
	346	int uneeded,
	347	int vtneeded,
	348	int additionalmemory,
	349	ref double[] w,
	350	ref double[,] u,
	351	ref double[,] vt)
[2154]	352	{
[2430]	353	bool result = new bool();
	354	double[] tauq = new double[0];
	355	double[] taup = new double[0];
	356	double[] tau = new double[0];
	357	double[] e = new double[0];
	358	double[] work = new double[0];
	359	double[,] t2 = new double[0,0];
	360	bool isupper = new bool();
	361	int minmn = 0;
	362	int ncu = 0;
	363	int nrvt = 0;
	364	int nru = 0;
	365	int ncvt = 0;
	366	int i = 0;
	367	int j = 0;
	368	int im1 = 0;
	369	double sm = 0;
	370
	371	a = (double[,])a.Clone();
	372
	373	result = true;
	374	if( m==0 \| n==0 )
	375	{
	376	return result;
	377	}
	378	System.Diagnostics.Debug.Assert(uneeded>=0 & uneeded<=2, "SVDDecomposition: wrong parameters!");
	379	System.Diagnostics.Debug.Assert(vtneeded>=0 & vtneeded<=2, "SVDDecomposition: wrong parameters!");
	380	System.Diagnostics.Debug.Assert(additionalmemory>=0 & additionalmemory<=2, "SVDDecomposition: wrong parameters!");
[2154]	381
	382	//
[2430]	383	// initialize
[2154]	384	//
[2430]	385	minmn = Math.Min(m, n);
	386	w = new double[minmn+1];
	387	ncu = 0;
	388	nru = 0;
	389	if( uneeded==1 )
	390	{
	391	nru = m;
	392	ncu = minmn;
	393	u = new double[nru+1, ncu+1];
	394	}
	395	if( uneeded==2 )
	396	{
	397	nru = m;
	398	ncu = m;
	399	u = new double[nru+1, ncu+1];
	400	}
	401	nrvt = 0;
	402	ncvt = 0;
	403	if( vtneeded==1 )
	404	{
	405	nrvt = minmn;
	406	ncvt = n;
	407	vt = new double[nrvt+1, ncvt+1];
	408	}
	409	if( vtneeded==2 )
	410	{
	411	nrvt = n;
	412	ncvt = n;
	413	vt = new double[nrvt+1, ncvt+1];
	414	}
[2154]	415
	416	//
[2430]	417	// M much larger than N
	418	// Use bidiagonal reduction with QR-decomposition
[2154]	419	//
[2563]	420	if( (double)(m)>(double)(1.6*n) )
[2154]	421	{
[2430]	422	if( uneeded==0 )
[2154]	423	{
[2430]	424
	425	//
	426	// No left singular vectors to be computed
	427	//
	428	qr.qrdecomposition(ref a, m, n, ref tau);
	429	for(i=2; i<=n; i++)
[2154]	430	{
[2430]	431	for(j=1; j<=i-1; j++)
	432	{
	433	a[i,j] = 0;
	434	}
[2154]	435	}
[2430]	436	bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
	437	bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
	438	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
	439	result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref a, 0, ref vt, ncvt);
	440	return result;
[2154]	441	}
[2430]	442	else
[2154]	443	{
[2430]	444
	445	//
	446	// Left singular vectors (may be full matrix U) to be computed
	447	//
	448	qr.qrdecomposition(ref a, m, n, ref tau);
	449	qr.unpackqfromqr(ref a, m, n, ref tau, ncu, ref u);
	450	for(i=2; i<=n; i++)
[2154]	451	{
[2430]	452	for(j=1; j<=i-1; j++)
	453	{
	454	a[i,j] = 0;
	455	}
[2154]	456	}
[2430]	457	bidiagonal.tobidiagonal(ref a, n, n, ref tauq, ref taup);
	458	bidiagonal.unpackptfrombidiagonal(ref a, n, n, ref taup, nrvt, ref vt);
	459	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, n, n, ref isupper, ref w, ref e);
	460	if( additionalmemory<1 )
	461	{
	462
	463	//
	464	// No additional memory can be used
	465	//
	466	bidiagonal.multiplybyqfrombidiagonal(ref a, n, n, ref tauq, ref u, m, n, true, false);
	467	result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, m, ref a, 0, ref vt, ncvt);
	468	}
	469	else
	470	{
	471
	472	//
	473	// Large U. Transforming intermediate matrix T2
	474	//
	475	work = new double[Math.Max(m, n)+1];
	476	bidiagonal.unpackqfrombidiagonal(ref a, n, n, ref tauq, n, ref t2);
	477	blas.copymatrix(ref u, 1, m, 1, n, ref a, 1, m, 1, n);
	478	blas.inplacetranspose(ref t2, 1, n, 1, n, ref work);
	479	result = bdsvd.bidiagonalsvddecomposition(ref w, e, n, isupper, false, ref u, 0, ref t2, n, ref vt, ncvt);
	480	blas.matrixmatrixmultiply(ref a, 1, m, 1, n, false, ref t2, 1, n, 1, n, true, 1.0, ref u, 1, m, 1, n, 0.0, ref work);
	481	}
	482	return result;
[2154]	483	}
[2430]	484	}
	485
	486	//
	487	// N much larger than M
	488	// Use bidiagonal reduction with LQ-decomposition
	489	//
[2563]	490	if( (double)(n)>(double)(1.6*m) )
[2430]	491	{
	492	if( vtneeded==0 )
[2154]	493	{
	494
	495	//
[2430]	496	// No right singular vectors to be computed
[2154]	497	//
[2430]	498	lq.lqdecomposition(ref a, m, n, ref tau);
	499	for(i=1; i<=m-1; i++)
	500	{
	501	for(j=i+1; j<=m; j++)
	502	{
	503	a[i,j] = 0;
	504	}
	505	}
	506	bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
	507	bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
	508	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
	509	work = new double[m+1];
	510	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
	511	result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, 0);
	512	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
	513	return result;
[2154]	514	}
	515	else
	516	{
	517
	518	//
[2430]	519	// Right singular vectors (may be full matrix VT) to be computed
[2154]	520	//
[2430]	521	lq.lqdecomposition(ref a, m, n, ref tau);
	522	lq.unpackqfromlq(ref a, m, n, ref tau, nrvt, ref vt);
	523	for(i=1; i<=m-1; i++)
	524	{
	525	for(j=i+1; j<=m; j++)
	526	{
	527	a[i,j] = 0;
	528	}
	529	}
	530	bidiagonal.tobidiagonal(ref a, m, m, ref tauq, ref taup);
	531	bidiagonal.unpackqfrombidiagonal(ref a, m, m, ref tauq, ncu, ref u);
	532	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, m, ref isupper, ref w, ref e);
[2154]	533	work = new double[Math.Max(m, n)+1];
[2430]	534	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
	535	if( additionalmemory<1 )
	536	{
	537
	538	//
	539	// No additional memory available
	540	//
	541	bidiagonal.multiplybypfrombidiagonal(ref a, m, m, ref taup, ref vt, m, n, false, true);
	542	result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref vt, n);
	543	}
	544	else
	545	{
	546
	547	//
	548	// Large VT. Transforming intermediate matrix T2
	549	//
	550	bidiagonal.unpackptfrombidiagonal(ref a, m, m, ref taup, m, ref t2);
	551	result = bdsvd.bidiagonalsvddecomposition(ref w, e, m, isupper, false, ref a, 0, ref u, nru, ref t2, m);
	552	blas.copymatrix(ref vt, 1, m, 1, n, ref a, 1, m, 1, n);
	553	blas.matrixmatrixmultiply(ref t2, 1, m, 1, m, false, ref a, 1, m, 1, n, false, 1.0, ref vt, 1, m, 1, n, 0.0, ref work);
	554	}
	555	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
	556	return result;
[2154]	557	}
	558	}
[2430]	559
	560	//
	561	// M<=N
	562	// We can use inplace transposition of U to get rid of columnwise operations
	563	//
	564	if( m<=n )
[2154]	565	{
[2430]	566	bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
	567	bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
	568	bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
	569	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
[2154]	570	work = new double[m+1];
	571	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
[2430]	572	result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref a, 0, ref u, nru, ref vt, ncvt);
[2154]	573	blas.inplacetranspose(ref u, 1, nru, 1, ncu, ref work);
	574	return result;
	575	}
[2430]	576
	577	//
	578	// Simple bidiagonal reduction
	579	//
	580	bidiagonal.tobidiagonal(ref a, m, n, ref tauq, ref taup);
	581	bidiagonal.unpackqfrombidiagonal(ref a, m, n, ref tauq, ncu, ref u);
	582	bidiagonal.unpackptfrombidiagonal(ref a, m, n, ref taup, nrvt, ref vt);
	583	bidiagonal.unpackdiagonalsfrombidiagonal(ref a, m, n, ref isupper, ref w, ref e);
	584	if( additionalmemory<2 \| uneeded==0 )
	585	{
	586
	587	//
	588	// We cant use additional memory or there is no need in such operations
	589	//
	590	result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, nru, ref a, 0, ref vt, ncvt);
	591	}
[2154]	592	else
	593	{
	594
	595	//
[2430]	596	// We can use additional memory
[2154]	597	//
[2430]	598	t2 = new double[minmn+1, m+1];
	599	blas.copyandtranspose(ref u, 1, m, 1, minmn, ref t2, 1, minmn, 1, m);
	600	result = bdsvd.bidiagonalsvddecomposition(ref w, e, minmn, isupper, false, ref u, 0, ref t2, m, ref vt, ncvt);
	601	blas.copyandtranspose(ref t2, 1, minmn, 1, m, ref u, 1, m, 1, minmn);
[2154]	602	}
	603	return result;
	604	}
	605	}
	606	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences