Context Navigation

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/bdsvd.cs @ 5955

Visit:

Last change on this file since 5955 was 2806, checked in by gkronber, 15 years ago
Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 50.6 KB

Line
1	/*************************************************************************
2	Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
3
4	Contributors:
5	* Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6	pseudocode.
7
8	See subroutines comments for additional copyrights.
9
10	>>> SOURCE LICENSE >>>
11	This program is free software; you can redistribute it and/or modify
12	it under the terms of the GNU General Public License as published by
13	the Free Software Foundation (www.fsf.org); either version 2 of the
14	License, or (at your option) any later version.
15
16	This program is distributed in the hope that it will be useful,
17	but WITHOUT ANY WARRANTY; without even the implied warranty of
18	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19	GNU General Public License for more details.
20
21	A copy of the GNU General Public License is available at
22	http://www.fsf.org/licensing/licenses
23
24	>>> END OF LICENSE >>>
25	*************************************************************************/
26
27	using System;
28
29	namespace alglib
30	{
31	public class bdsvd
32	{
33	/*************************************************************************
34	Singular value decomposition of a bidiagonal matrix (extended algorithm)
35
36	The algorithm performs the singular value decomposition of a bidiagonal
37	matrix B (upper or lower) representing it as B = QSP^T, where Q and P -
38	orthogonal matrices, S - diagonal matrix with non-negative elements on the
39	main diagonal, in descending order.
40
41	The algorithm finds singular values. In addition, the algorithm can
42	calculate matrices Q and P (more precisely, not the matrices, but their
43	product with given matrices U and VT - UQ and (P^T)VT)). Of course,
44	matrices U and VT can be of any type, including identity. Furthermore, the
45	algorithm can calculate Q'*C (this product is calculated more effectively
46	than U*Q, because this calculation operates with rows instead of matrix
47	columns).
48
49	The feature of the algorithm is its ability to find all singular values
50	including those which are arbitrarily close to 0 with relative accuracy
51	close to machine precision. If the parameter IsFractionalAccuracyRequired
52	is set to True, all singular values will have high relative accuracy close
53	to machine precision. If the parameter is set to False, only the biggest
54	singular value will have relative accuracy close to machine precision.
55	The absolute error of other singular values is equal to the absolute error
56	of the biggest singular value.
57
58	Input parameters:
59	D - main diagonal of matrix B.
60	Array whose index ranges within [0..N-1].
61	E - superdiagonal (or subdiagonal) of matrix B.
62	Array whose index ranges within [0..N-2].
63	N - size of matrix B.
64	IsUpper - True, if the matrix is upper bidiagonal.
65	IsFractionalAccuracyRequired -
66	accuracy to search singular values with.
67	U - matrix to be multiplied by Q.
68	Array whose indexes range within [0..NRU-1, 0..N-1].
69	The matrix can be bigger, in that case only the submatrix
70	[0..NRU-1, 0..N-1] will be multiplied by Q.
71	NRU - number of rows in matrix U.
72	C - matrix to be multiplied by Q'.
73	Array whose indexes range within [0..N-1, 0..NCC-1].
74	The matrix can be bigger, in that case only the submatrix
75	[0..N-1, 0..NCC-1] will be multiplied by Q'.
76	NCC - number of columns in matrix C.
77	VT - matrix to be multiplied by P^T.
78	Array whose indexes range within [0..N-1, 0..NCVT-1].
79	The matrix can be bigger, in that case only the submatrix
80	[0..N-1, 0..NCVT-1] will be multiplied by P^T.
81	NCVT - number of columns in matrix VT.
82
83	Output parameters:
84	D - singular values of matrix B in descending order.
85	U - if NRU>0, contains matrix U*Q.
86	VT - if NCVT>0, contains matrix (P^T)*VT.
87	C - if NCC>0, contains matrix Q'*C.
88
89	Result:
90	True, if the algorithm has converged.
91	False, if the algorithm hasn't converged (rare case).
92
93	Additional information:
94	The type of convergence is controlled by the internal parameter TOL.
95	If the parameter is greater than 0, the singular values will have
96	relative accuracy TOL. If TOL<0, the singular values will have
97	absolute accuracy ABS(TOL)*norm(B).
98	By default, \|TOL\| falls within the range of 10Epsilon and 100Epsilon,
99	where Epsilon is the machine precision. It is not recommended to use
100	TOL less than 10*Epsilon since this will considerably slow down the
101	algorithm and may not lead to error decreasing.
102	History:
103	* 31 March, 2007.
104	changed MAXITR from 6 to 12.
105
106	-- LAPACK routine (version 3.0) --
107	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
108	Courant Institute, Argonne National Lab, and Rice University
109	October 31, 1999.
110	*************************************************************************/
111	public static bool rmatrixbdsvd(ref double[] d,
112	double[] e,
113	int n,
114	bool isupper,
115	bool isfractionalaccuracyrequired,
116	ref double[,] u,
117	int nru,
118	ref double[,] c,
119	int ncc,
120	ref double[,] vt,
121	int ncvt)
122	{
123	bool result = new bool();
124	double[] d1 = new double[0];
125	double[] e1 = new double[0];
126	int i_ = 0;
127	int i1_ = 0;
128
129	e = (double[])e.Clone();
130
131	d1 = new double[n+1];
132	i1_ = (0) - (1);
133	for(i_=1; i_<=n;i_++)
134	{
135	d1[i_] = d[i_+i1_];
136	}
137	if( n>1 )
138	{
139	e1 = new double[n-1+1];
140	i1_ = (0) - (1);
141	for(i_=1; i_<=n-1;i_++)
142	{
143	e1[i_] = e[i_+i1_];
144	}
145	}
146	result = bidiagonalsvddecompositioninternal(ref d1, e1, n, isupper, isfractionalaccuracyrequired, ref u, 0, nru, ref c, 0, ncc, ref vt, 0, ncvt);
147	i1_ = (1) - (0);
148	for(i_=0; i_<=n-1;i_++)
149	{
150	d[i_] = d1[i_+i1_];
151	}
152	return result;
153	}
154
155
156	public static bool bidiagonalsvddecomposition(ref double[] d,
157	double[] e,
158	int n,
159	bool isupper,
160	bool isfractionalaccuracyrequired,
161	ref double[,] u,
162	int nru,
163	ref double[,] c,
164	int ncc,
165	ref double[,] vt,
166	int ncvt)
167	{
168	bool result = new bool();
169
170	e = (double[])e.Clone();
171
172	result = bidiagonalsvddecompositioninternal(ref d, e, n, isupper, isfractionalaccuracyrequired, ref u, 1, nru, ref c, 1, ncc, ref vt, 1, ncvt);
173	return result;
174	}
175
176
177	/*************************************************************************
178	Internal working subroutine for bidiagonal decomposition
179	*************************************************************************/
180	private static bool bidiagonalsvddecompositioninternal(ref double[] d,
181	double[] e,
182	int n,
183	bool isupper,
184	bool isfractionalaccuracyrequired,
185	ref double[,] u,
186	int ustart,
187	int nru,
188	ref double[,] c,
189	int cstart,
190	int ncc,
191	ref double[,] vt,
192	int vstart,
193	int ncvt)
194	{
195	bool result = new bool();
196	int i = 0;
197	int idir = 0;
198	int isub = 0;
199	int iter = 0;
200	int j = 0;
201	int ll = 0;
202	int lll = 0;
203	int m = 0;
204	int maxit = 0;
205	int oldll = 0;
206	int oldm = 0;
207	double abse = 0;
208	double abss = 0;
209	double cosl = 0;
210	double cosr = 0;
211	double cs = 0;
212	double eps = 0;
213	double f = 0;
214	double g = 0;
215	double h = 0;
216	double mu = 0;
217	double oldcs = 0;
218	double oldsn = 0;
219	double r = 0;
220	double shift = 0;
221	double sigmn = 0;
222	double sigmx = 0;
223	double sinl = 0;
224	double sinr = 0;
225	double sll = 0;
226	double smax = 0;
227	double smin = 0;
228	double sminl = 0;
229	double sminlo = 0;
230	double sminoa = 0;
231	double sn = 0;
232	double thresh = 0;
233	double tol = 0;
234	double tolmul = 0;
235	double unfl = 0;
236	double[] work0 = new double[0];
237	double[] work1 = new double[0];
238	double[] work2 = new double[0];
239	double[] work3 = new double[0];
240	int maxitr = 0;
241	bool matrixsplitflag = new bool();
242	bool iterflag = new bool();
243	double[] utemp = new double[0];
244	double[] vttemp = new double[0];
245	double[] ctemp = new double[0];
246	double[] etemp = new double[0];
247	bool rightside = new bool();
248	bool fwddir = new bool();
249	double tmp = 0;
250	int mm1 = 0;
251	int mm0 = 0;
252	bool bchangedir = new bool();
253	int uend = 0;
254	int cend = 0;
255	int vend = 0;
256	int i_ = 0;
257
258	e = (double[])e.Clone();
259
260	result = true;
261	if( n==0 )
262	{
263	return result;
264	}
265	if( n==1 )
266	{
267	if( (double)(d[1])<(double)(0) )
268	{
269	d[1] = -d[1];
270	if( ncvt>0 )
271	{
272	for(i_=vstart; i_<=vstart+ncvt-1;i_++)
273	{
274	vt[vstart,i_] = -1*vt[vstart,i_];
275	}
276	}
277	}
278	return result;
279	}
280
281	//
282	// init
283	//
284	work0 = new double[n-1+1];
285	work1 = new double[n-1+1];
286	work2 = new double[n-1+1];
287	work3 = new double[n-1+1];
288	uend = ustart+Math.Max(nru-1, 0);
289	vend = vstart+Math.Max(ncvt-1, 0);
290	cend = cstart+Math.Max(ncc-1, 0);
291	utemp = new double[uend+1];
292	vttemp = new double[vend+1];
293	ctemp = new double[cend+1];
294	maxitr = 12;
295	rightside = true;
296	fwddir = true;
297
298	//
299	// resize E from N-1 to N
300	//
301	etemp = new double[n+1];
302	for(i=1; i<=n-1; i++)
303	{
304	etemp[i] = e[i];
305	}
306	e = new double[n+1];
307	for(i=1; i<=n-1; i++)
308	{
309	e[i] = etemp[i];
310	}
311	e[n] = 0;
312	idir = 0;
313
314	//
315	// Get machine constants
316	//
317	eps = AP.Math.MachineEpsilon;
318	unfl = AP.Math.MinRealNumber;
319
320	//
321	// If matrix lower bidiagonal, rotate to be upper bidiagonal
322	// by applying Givens rotations on the left
323	//
324	if( !isupper )
325	{
326	for(i=1; i<=n-1; i++)
327	{
328	rotations.generaterotation(d[i], e[i], ref cs, ref sn, ref r);
329	d[i] = r;
330	e[i] = sn*d[i+1];
331	d[i+1] = cs*d[i+1];
332	work0[i] = cs;
333	work1[i] = sn;
334	}
335
336	//
337	// Update singular vectors if desired
338	//
339	if( nru>0 )
340	{
341	rotations.applyrotationsfromtheright(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, ref work0, ref work1, ref u, ref utemp);
342	}
343	if( ncc>0 )
344	{
345	rotations.applyrotationsfromtheleft(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
346	}
347	}
348
349	//
350	// Compute singular values to relative accuracy TOL
351	// (By setting TOL to be negative, algorithm will compute
352	// singular values to absolute accuracy ABS(TOL)*norm(input matrix))
353	//
354	tolmul = Math.Max(10, Math.Min(100, Math.Pow(eps, -0.125)));
355	tol = tolmul*eps;
356	if( !isfractionalaccuracyrequired )
357	{
358	tol = -tol;
359	}
360
361	//
362	// Compute approximate maximum, minimum singular values
363	//
364	smax = 0;
365	for(i=1; i<=n; i++)
366	{
367	smax = Math.Max(smax, Math.Abs(d[i]));
368	}
369	for(i=1; i<=n-1; i++)
370	{
371	smax = Math.Max(smax, Math.Abs(e[i]));
372	}
373	sminl = 0;
374	if( (double)(tol)>=(double)(0) )
375	{
376
377	//
378	// Relative accuracy desired
379	//
380	sminoa = Math.Abs(d[1]);
381	if( (double)(sminoa)!=(double)(0) )
382	{
383	mu = sminoa;
384	for(i=2; i<=n; i++)
385	{
386	mu = Math.Abs(d[i])*(mu/(mu+Math.Abs(e[i-1])));
387	sminoa = Math.Min(sminoa, mu);
388	if( (double)(sminoa)==(double)(0) )
389	{
390	break;
391	}
392	}
393	}
394	sminoa = sminoa/Math.Sqrt(n);
395	thresh = Math.Max(tolsminoa, maxitrnnunfl);
396	}
397	else
398	{
399
400	//
401	// Absolute accuracy desired
402	//
403	thresh = Math.Max(Math.Abs(tol)smax, maxitrnnunfl);
404	}
405
406	//
407	// Prepare for main iteration loop for the singular values
408	// (MAXIT is the maximum number of passes through the inner
409	// loop permitted before nonconvergence signalled.)
410	//
411	maxit = maxitrnn;
412	iter = 0;
413	oldll = -1;
414	oldm = -1;
415
416	//
417	// M points to last element of unconverged part of matrix
418	//
419	m = n;
420
421	//
422	// Begin main iteration loop
423	//
424	while( true )
425	{
426
427	//
428	// Check for convergence or exceeding iteration count
429	//
430	if( m<=1 )
431	{
432	break;
433	}
434	if( iter>maxit )
435	{
436	result = false;
437	return result;
438	}
439
440	//
441	// Find diagonal block of matrix to work on
442	//
443	if( (double)(tol)<(double)(0) & (double)(Math.Abs(d[m]))<=(double)(thresh) )
444	{
445	d[m] = 0;
446	}
447	smax = Math.Abs(d[m]);
448	smin = smax;
449	matrixsplitflag = false;
450	for(lll=1; lll<=m-1; lll++)
451	{
452	ll = m-lll;
453	abss = Math.Abs(d[ll]);
454	abse = Math.Abs(e[ll]);
455	if( (double)(tol)<(double)(0) & (double)(abss)<=(double)(thresh) )
456	{
457	d[ll] = 0;
458	}
459	if( (double)(abse)<=(double)(thresh) )
460	{
461	matrixsplitflag = true;
462	break;
463	}
464	smin = Math.Min(smin, abss);
465	smax = Math.Max(smax, Math.Max(abss, abse));
466	}
467	if( !matrixsplitflag )
468	{
469	ll = 0;
470	}
471	else
472	{
473
474	//
475	// Matrix splits since E(LL) = 0
476	//
477	e[ll] = 0;
478	if( ll==m-1 )
479	{
480
481	//
482	// Convergence of bottom singular value, return to top of loop
483	//
484	m = m-1;
485	continue;
486	}
487	}
488	ll = ll+1;
489
490	//
491	// E(LL) through E(M-1) are nonzero, E(LL-1) is zero
492	//
493	if( ll==m-1 )
494	{
495
496	//
497	// 2 by 2 block, handle separately
498	//
499	svdv2x2(d[m-1], e[m-1], d[m], ref sigmn, ref sigmx, ref sinr, ref cosr, ref sinl, ref cosl);
500	d[m-1] = sigmx;
501	e[m-1] = 0;
502	d[m] = sigmn;
503
504	//
505	// Compute singular vectors, if desired
506	//
507	if( ncvt>0 )
508	{
509	mm0 = m+(vstart-1);
510	mm1 = m-1+(vstart-1);
511	for(i_=vstart; i_<=vend;i_++)
512	{
513	vttemp[i_] = cosr*vt[mm1,i_];
514	}
515	for(i_=vstart; i_<=vend;i_++)
516	{
517	vttemp[i_] = vttemp[i_] + sinr*vt[mm0,i_];
518	}
519	for(i_=vstart; i_<=vend;i_++)
520	{
521	vt[mm0,i_] = cosr*vt[mm0,i_];
522	}
523	for(i_=vstart; i_<=vend;i_++)
524	{
525	vt[mm0,i_] = vt[mm0,i_] - sinr*vt[mm1,i_];
526	}
527	for(i_=vstart; i_<=vend;i_++)
528	{
529	vt[mm1,i_] = vttemp[i_];
530	}
531	}
532	if( nru>0 )
533	{
534	mm0 = m+ustart-1;
535	mm1 = m-1+ustart-1;
536	for(i_=ustart; i_<=uend;i_++)
537	{
538	utemp[i_] = cosl*u[i_,mm1];
539	}
540	for(i_=ustart; i_<=uend;i_++)
541	{
542	utemp[i_] = utemp[i_] + sinl*u[i_,mm0];
543	}
544	for(i_=ustart; i_<=uend;i_++)
545	{
546	u[i_,mm0] = cosl*u[i_,mm0];
547	}
548	for(i_=ustart; i_<=uend;i_++)
549	{
550	u[i_,mm0] = u[i_,mm0] - sinl*u[i_,mm1];
551	}
552	for(i_=ustart; i_<=uend;i_++)
553	{
554	u[i_,mm1] = utemp[i_];
555	}
556	}
557	if( ncc>0 )
558	{
559	mm0 = m+cstart-1;
560	mm1 = m-1+cstart-1;
561	for(i_=cstart; i_<=cend;i_++)
562	{
563	ctemp[i_] = cosl*c[mm1,i_];
564	}
565	for(i_=cstart; i_<=cend;i_++)
566	{
567	ctemp[i_] = ctemp[i_] + sinl*c[mm0,i_];
568	}
569	for(i_=cstart; i_<=cend;i_++)
570	{
571	c[mm0,i_] = cosl*c[mm0,i_];
572	}
573	for(i_=cstart; i_<=cend;i_++)
574	{
575	c[mm0,i_] = c[mm0,i_] - sinl*c[mm1,i_];
576	}
577	for(i_=cstart; i_<=cend;i_++)
578	{
579	c[mm1,i_] = ctemp[i_];
580	}
581	}
582	m = m-2;
583	continue;
584	}
585
586	//
587	// If working on new submatrix, choose shift direction
588	// (from larger end diagonal element towards smaller)
589	//
590	// Previously was
591	// "if (LL>OLDM) or (M<OLDLL) then"
592	// fixed thanks to Michael Rolle < m@rolle.name >
593	// Very strange that LAPACK still contains it.
594	//
595	bchangedir = false;
596	if( idir==1 & (double)(Math.Abs(d[ll]))<(double)(1.0E-3*Math.Abs(d[m])) )
597	{
598	bchangedir = true;
599	}
600	if( idir==2 & (double)(Math.Abs(d[m]))<(double)(1.0E-3*Math.Abs(d[ll])) )
601	{
602	bchangedir = true;
603	}
604	if( ll!=oldll \| m!=oldm \| bchangedir )
605	{
606	if( (double)(Math.Abs(d[ll]))>=(double)(Math.Abs(d[m])) )
607	{
608
609	//
610	// Chase bulge from top (big end) to bottom (small end)
611	//
612	idir = 1;
613	}
614	else
615	{
616
617	//
618	// Chase bulge from bottom (big end) to top (small end)
619	//
620	idir = 2;
621	}
622	}
623
624	//
625	// Apply convergence tests
626	//
627	if( idir==1 )
628	{
629
630	//
631	// Run convergence test in forward direction
632	// First apply standard test to bottom of matrix
633	//
634	if( (double)(Math.Abs(e[m-1]))<=(double)(Math.Abs(tol)*Math.Abs(d[m])) \| (double)(tol)<(double)(0) & (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
635	{
636	e[m-1] = 0;
637	continue;
638	}
639	if( (double)(tol)>=(double)(0) )
640	{
641
642	//
643	// If relative accuracy desired,
644	// apply convergence criterion forward
645	//
646	mu = Math.Abs(d[ll]);
647	sminl = mu;
648	iterflag = false;
649	for(lll=ll; lll<=m-1; lll++)
650	{
651	if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
652	{
653	e[lll] = 0;
654	iterflag = true;
655	break;
656	}
657	sminlo = sminl;
658	mu = Math.Abs(d[lll+1])*(mu/(mu+Math.Abs(e[lll])));
659	sminl = Math.Min(sminl, mu);
660	}
661	if( iterflag )
662	{
663	continue;
664	}
665	}
666	}
667	else
668	{
669
670	//
671	// Run convergence test in backward direction
672	// First apply standard test to top of matrix
673	//
674	if( (double)(Math.Abs(e[ll]))<=(double)(Math.Abs(tol)*Math.Abs(d[ll])) \| (double)(tol)<(double)(0) & (double)(Math.Abs(e[ll]))<=(double)(thresh) )
675	{
676	e[ll] = 0;
677	continue;
678	}
679	if( (double)(tol)>=(double)(0) )
680	{
681
682	//
683	// If relative accuracy desired,
684	// apply convergence criterion backward
685	//
686	mu = Math.Abs(d[m]);
687	sminl = mu;
688	iterflag = false;
689	for(lll=m-1; lll>=ll; lll--)
690	{
691	if( (double)(Math.Abs(e[lll]))<=(double)(tol*mu) )
692	{
693	e[lll] = 0;
694	iterflag = true;
695	break;
696	}
697	sminlo = sminl;
698	mu = Math.Abs(d[lll])*(mu/(mu+Math.Abs(e[lll])));
699	sminl = Math.Min(sminl, mu);
700	}
701	if( iterflag )
702	{
703	continue;
704	}
705	}
706	}
707	oldll = ll;
708	oldm = m;
709
710	//
711	// Compute shift. First, test if shifting would ruin relative
712	// accuracy, and if so set the shift to zero.
713	//
714	if( (double)(tol)>=(double)(0) & (double)(ntol(sminl/smax))<=(double)(Math.Max(eps, 0.01*tol)) )
715	{
716
717	//
718	// Use a zero shift to avoid loss of relative accuracy
719	//
720	shift = 0;
721	}
722	else
723	{
724
725	//
726	// Compute the shift from 2-by-2 block at end of matrix
727	//
728	if( idir==1 )
729	{
730	sll = Math.Abs(d[ll]);
731	svd2x2(d[m-1], e[m-1], d[m], ref shift, ref r);
732	}
733	else
734	{
735	sll = Math.Abs(d[m]);
736	svd2x2(d[ll], e[ll], d[ll+1], ref shift, ref r);
737	}
738
739	//
740	// Test if shift negligible, and if so set to zero
741	//
742	if( (double)(sll)>(double)(0) )
743	{
744	if( (double)(AP.Math.Sqr(shift/sll))<(double)(eps) )
745	{
746	shift = 0;
747	}
748	}
749	}
750
751	//
752	// Increment iteration count
753	//
754	iter = iter+m-ll;
755
756	//
757	// If SHIFT = 0, do simplified QR iteration
758	//
759	if( (double)(shift)==(double)(0) )
760	{
761	if( idir==1 )
762	{
763
764	//
765	// Chase bulge from top to bottom
766	// Save cosines and sines for later singular vector updates
767	//
768	cs = 1;
769	oldcs = 1;
770	for(i=ll; i<=m-1; i++)
771	{
772	rotations.generaterotation(d[i]*cs, e[i], ref cs, ref sn, ref r);
773	if( i>ll )
774	{
775	e[i-1] = oldsn*r;
776	}
777	rotations.generaterotation(oldcsr, d[i+1]sn, ref oldcs, ref oldsn, ref tmp);
778	d[i] = tmp;
779	work0[i-ll+1] = cs;
780	work1[i-ll+1] = sn;
781	work2[i-ll+1] = oldcs;
782	work3[i-ll+1] = oldsn;
783	}
784	h = d[m]*cs;
785	d[m] = h*oldcs;
786	e[m-1] = h*oldsn;
787
788	//
789	// Update singular vectors
790	//
791	if( ncvt>0 )
792	{
793	rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
794	}
795	if( nru>0 )
796	{
797	rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
798	}
799	if( ncc>0 )
800	{
801	rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
802	}
803
804	//
805	// Test convergence
806	//
807	if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
808	{
809	e[m-1] = 0;
810	}
811	}
812	else
813	{
814
815	//
816	// Chase bulge from bottom to top
817	// Save cosines and sines for later singular vector updates
818	//
819	cs = 1;
820	oldcs = 1;
821	for(i=m; i>=ll+1; i--)
822	{
823	rotations.generaterotation(d[i]*cs, e[i-1], ref cs, ref sn, ref r);
824	if( i<m )
825	{
826	e[i] = oldsn*r;
827	}
828	rotations.generaterotation(oldcsr, d[i-1]sn, ref oldcs, ref oldsn, ref tmp);
829	d[i] = tmp;
830	work0[i-ll] = cs;
831	work1[i-ll] = -sn;
832	work2[i-ll] = oldcs;
833	work3[i-ll] = -oldsn;
834	}
835	h = d[ll]*cs;
836	d[ll] = h*oldcs;
837	e[ll] = h*oldsn;
838
839	//
840	// Update singular vectors
841	//
842	if( ncvt>0 )
843	{
844	rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
845	}
846	if( nru>0 )
847	{
848	rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
849	}
850	if( ncc>0 )
851	{
852	rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
853	}
854
855	//
856	// Test convergence
857	//
858	if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
859	{
860	e[ll] = 0;
861	}
862	}
863	}
864	else
865	{
866
867	//
868	// Use nonzero shift
869	//
870	if( idir==1 )
871	{
872
873	//
874	// Chase bulge from top to bottom
875	// Save cosines and sines for later singular vector updates
876	//
877	f = (Math.Abs(d[ll])-shift)*(extsignbdsqr(1, d[ll])+shift/d[ll]);
878	g = e[ll];
879	for(i=ll; i<=m-1; i++)
880	{
881	rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
882	if( i>ll )
883	{
884	e[i-1] = r;
885	}
886	f = cosrd[i]+sinre[i];
887	e[i] = cosre[i]-sinrd[i];
888	g = sinr*d[i+1];
889	d[i+1] = cosr*d[i+1];
890	rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
891	d[i] = r;
892	f = cosle[i]+sinld[i+1];
893	d[i+1] = cosld[i+1]-sinle[i];
894	if( i<m-1 )
895	{
896	g = sinl*e[i+1];
897	e[i+1] = cosl*e[i+1];
898	}
899	work0[i-ll+1] = cosr;
900	work1[i-ll+1] = sinr;
901	work2[i-ll+1] = cosl;
902	work3[i-ll+1] = sinl;
903	}
904	e[m-1] = f;
905
906	//
907	// Update singular vectors
908	//
909	if( ncvt>0 )
910	{
911	rotations.applyrotationsfromtheleft(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work0, ref work1, ref vt, ref vttemp);
912	}
913	if( nru>0 )
914	{
915	rotations.applyrotationsfromtheright(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work2, ref work3, ref u, ref utemp);
916	}
917	if( ncc>0 )
918	{
919	rotations.applyrotationsfromtheleft(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work2, ref work3, ref c, ref ctemp);
920	}
921
922	//
923	// Test convergence
924	//
925	if( (double)(Math.Abs(e[m-1]))<=(double)(thresh) )
926	{
927	e[m-1] = 0;
928	}
929	}
930	else
931	{
932
933	//
934	// Chase bulge from bottom to top
935	// Save cosines and sines for later singular vector updates
936	//
937	f = (Math.Abs(d[m])-shift)*(extsignbdsqr(1, d[m])+shift/d[m]);
938	g = e[m-1];
939	for(i=m; i>=ll+1; i--)
940	{
941	rotations.generaterotation(f, g, ref cosr, ref sinr, ref r);
942	if( i<m )
943	{
944	e[i] = r;
945	}
946	f = cosrd[i]+sinre[i-1];
947	e[i-1] = cosre[i-1]-sinrd[i];
948	g = sinr*d[i-1];
949	d[i-1] = cosr*d[i-1];
950	rotations.generaterotation(f, g, ref cosl, ref sinl, ref r);
951	d[i] = r;
952	f = cosle[i-1]+sinld[i-1];
953	d[i-1] = cosld[i-1]-sinle[i-1];
954	if( i>ll+1 )
955	{
956	g = sinl*e[i-2];
957	e[i-2] = cosl*e[i-2];
958	}
959	work0[i-ll] = cosr;
960	work1[i-ll] = -sinr;
961	work2[i-ll] = cosl;
962	work3[i-ll] = -sinl;
963	}
964	e[ll] = f;
965
966	//
967	// Test convergence
968	//
969	if( (double)(Math.Abs(e[ll]))<=(double)(thresh) )
970	{
971	e[ll] = 0;
972	}
973
974	//
975	// Update singular vectors if desired
976	//
977	if( ncvt>0 )
978	{
979	rotations.applyrotationsfromtheleft(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, ref work2, ref work3, ref vt, ref vttemp);
980	}
981	if( nru>0 )
982	{
983	rotations.applyrotationsfromtheright(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, ref work0, ref work1, ref u, ref utemp);
984	}
985	if( ncc>0 )
986	{
987	rotations.applyrotationsfromtheleft(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, ref work0, ref work1, ref c, ref ctemp);
988	}
989	}
990	}
991
992	//
993	// QR iteration finished, go back and check convergence
994	//
995	continue;
996	}
997
998	//
999	// All singular values converged, so make them positive
1000	//
1001	for(i=1; i<=n; i++)
1002	{
1003	if( (double)(d[i])<(double)(0) )
1004	{
1005	d[i] = -d[i];
1006
1007	//
1008	// Change sign of singular vectors, if desired
1009	//
1010	if( ncvt>0 )
1011	{
1012	for(i_=vstart; i_<=vend;i_++)
1013	{
1014	vt[i+vstart-1,i_] = -1*vt[i+vstart-1,i_];
1015	}
1016	}
1017	}
1018	}
1019
1020	//
1021	// Sort the singular values into decreasing order (insertion sort on
1022	// singular values, but only one transposition per singular vector)
1023	//
1024	for(i=1; i<=n-1; i++)
1025	{
1026
1027	//
1028	// Scan for smallest D(I)
1029	//
1030	isub = 1;
1031	smin = d[1];
1032	for(j=2; j<=n+1-i; j++)
1033	{
1034	if( (double)(d[j])<=(double)(smin) )
1035	{
1036	isub = j;
1037	smin = d[j];
1038	}
1039	}
1040	if( isub!=n+1-i )
1041	{
1042
1043	//
1044	// Swap singular values and vectors
1045	//
1046	d[isub] = d[n+1-i];
1047	d[n+1-i] = smin;
1048	if( ncvt>0 )
1049	{
1050	j = n+1-i;
1051	for(i_=vstart; i_<=vend;i_++)
1052	{
1053	vttemp[i_] = vt[isub+vstart-1,i_];
1054	}
1055	for(i_=vstart; i_<=vend;i_++)
1056	{
1057	vt[isub+vstart-1,i_] = vt[j+vstart-1,i_];
1058	}
1059	for(i_=vstart; i_<=vend;i_++)
1060	{
1061	vt[j+vstart-1,i_] = vttemp[i_];
1062	}
1063	}
1064	if( nru>0 )
1065	{
1066	j = n+1-i;
1067	for(i_=ustart; i_<=uend;i_++)
1068	{
1069	utemp[i_] = u[i_,isub+ustart-1];
1070	}
1071	for(i_=ustart; i_<=uend;i_++)
1072	{
1073	u[i_,isub+ustart-1] = u[i_,j+ustart-1];
1074	}
1075	for(i_=ustart; i_<=uend;i_++)
1076	{
1077	u[i_,j+ustart-1] = utemp[i_];
1078	}
1079	}
1080	if( ncc>0 )
1081	{
1082	j = n+1-i;
1083	for(i_=cstart; i_<=cend;i_++)
1084	{
1085	ctemp[i_] = c[isub+cstart-1,i_];
1086	}
1087	for(i_=cstart; i_<=cend;i_++)
1088	{
1089	c[isub+cstart-1,i_] = c[j+cstart-1,i_];
1090	}
1091	for(i_=cstart; i_<=cend;i_++)
1092	{
1093	c[j+cstart-1,i_] = ctemp[i_];
1094	}
1095	}
1096	}
1097	}
1098	return result;
1099	}
1100
1101
1102	private static double extsignbdsqr(double a,
1103	double b)
1104	{
1105	double result = 0;
1106
1107	if( (double)(b)>=(double)(0) )
1108	{
1109	result = Math.Abs(a);
1110	}
1111	else
1112	{
1113	result = -Math.Abs(a);
1114	}
1115	return result;
1116	}
1117
1118
1119	private static void svd2x2(double f,
1120	double g,
1121	double h,
1122	ref double ssmin,
1123	ref double ssmax)
1124	{
1125	double aas = 0;
1126	double at = 0;
1127	double au = 0;
1128	double c = 0;
1129	double fa = 0;
1130	double fhmn = 0;
1131	double fhmx = 0;
1132	double ga = 0;
1133	double ha = 0;
1134
1135	fa = Math.Abs(f);
1136	ga = Math.Abs(g);
1137	ha = Math.Abs(h);
1138	fhmn = Math.Min(fa, ha);
1139	fhmx = Math.Max(fa, ha);
1140	if( (double)(fhmn)==(double)(0) )
1141	{
1142	ssmin = 0;
1143	if( (double)(fhmx)==(double)(0) )
1144	{
1145	ssmax = ga;
1146	}
1147	else
1148	{
1149	ssmax = Math.Max(fhmx, ga)*Math.Sqrt(1+AP.Math.Sqr(Math.Min(fhmx, ga)/Math.Max(fhmx, ga)));
1150	}
1151	}
1152	else
1153	{
1154	if( (double)(ga)<(double)(fhmx) )
1155	{
1156	aas = 1+fhmn/fhmx;
1157	at = (fhmx-fhmn)/fhmx;
1158	au = AP.Math.Sqr(ga/fhmx);
1159	c = 2/(Math.Sqrt(aasaas+au)+Math.Sqrt(atat+au));
1160	ssmin = fhmn*c;
1161	ssmax = fhmx/c;
1162	}
1163	else
1164	{
1165	au = fhmx/ga;
1166	if( (double)(au)==(double)(0) )
1167	{
1168
1169	//
1170	// Avoid possible harmful underflow if exponent range
1171	// asymmetric (true SSMIN may not underflow even if
1172	// AU underflows)
1173	//
1174	ssmin = fhmn*fhmx/ga;
1175	ssmax = ga;
1176	}
1177	else
1178	{
1179	aas = 1+fhmn/fhmx;
1180	at = (fhmx-fhmn)/fhmx;
1181	c = 1/(Math.Sqrt(1+AP.Math.Sqr(aasau))+Math.Sqrt(1+AP.Math.Sqr(atau)));
1182	ssmin = fhmncau;
1183	ssmin = ssmin+ssmin;
1184	ssmax = ga/(c+c);
1185	}
1186	}
1187	}
1188	}
1189
1190
1191	private static void svdv2x2(double f,
1192	double g,
1193	double h,
1194	ref double ssmin,
1195	ref double ssmax,
1196	ref double snr,
1197	ref double csr,
1198	ref double snl,
1199	ref double csl)
1200	{
1201	bool gasmal = new bool();
1202	bool swp = new bool();
1203	int pmax = 0;
1204	double a = 0;
1205	double clt = 0;
1206	double crt = 0;
1207	double d = 0;
1208	double fa = 0;
1209	double ft = 0;
1210	double ga = 0;
1211	double gt = 0;
1212	double ha = 0;
1213	double ht = 0;
1214	double l = 0;
1215	double m = 0;
1216	double mm = 0;
1217	double r = 0;
1218	double s = 0;
1219	double slt = 0;
1220	double srt = 0;
1221	double t = 0;
1222	double temp = 0;
1223	double tsign = 0;
1224	double tt = 0;
1225	double v = 0;
1226
1227	ft = f;
1228	fa = Math.Abs(ft);
1229	ht = h;
1230	ha = Math.Abs(h);
1231
1232	//
1233	// PMAX points to the maximum absolute element of matrix
1234	// PMAX = 1 if F largest in absolute values
1235	// PMAX = 2 if G largest in absolute values
1236	// PMAX = 3 if H largest in absolute values
1237	//
1238	pmax = 1;
1239	swp = (double)(ha)>(double)(fa);
1240	if( swp )
1241	{
1242
1243	//
1244	// Now FA .ge. HA
1245	//
1246	pmax = 3;
1247	temp = ft;
1248	ft = ht;
1249	ht = temp;
1250	temp = fa;
1251	fa = ha;
1252	ha = temp;
1253	}
1254	gt = g;
1255	ga = Math.Abs(gt);
1256	if( (double)(ga)==(double)(0) )
1257	{
1258
1259	//
1260	// Diagonal matrix
1261	//
1262	ssmin = ha;
1263	ssmax = fa;
1264	clt = 1;
1265	crt = 1;
1266	slt = 0;
1267	srt = 0;
1268	}
1269	else
1270	{
1271	gasmal = true;
1272	if( (double)(ga)>(double)(fa) )
1273	{
1274	pmax = 2;
1275	if( (double)(fa/ga)<(double)(AP.Math.MachineEpsilon) )
1276	{
1277
1278	//
1279	// Case of very large GA
1280	//
1281	gasmal = false;
1282	ssmax = ga;
1283	if( (double)(ha)>(double)(1) )
1284	{
1285	v = ga/ha;
1286	ssmin = fa/v;
1287	}
1288	else
1289	{
1290	v = fa/ga;
1291	ssmin = v*ha;
1292	}
1293	clt = 1;
1294	slt = ht/gt;
1295	srt = 1;
1296	crt = ft/gt;
1297	}
1298	}
1299	if( gasmal )
1300	{
1301
1302	//
1303	// Normal case
1304	//
1305	d = fa-ha;
1306	if( (double)(d)==(double)(fa) )
1307	{
1308	l = 1;
1309	}
1310	else
1311	{
1312	l = d/fa;
1313	}
1314	m = gt/ft;
1315	t = 2-l;
1316	mm = m*m;
1317	tt = t*t;
1318	s = Math.Sqrt(tt+mm);
1319	if( (double)(l)==(double)(0) )
1320	{
1321	r = Math.Abs(m);
1322	}
1323	else
1324	{
1325	r = Math.Sqrt(l*l+mm);
1326	}
1327	a = 0.5*(s+r);
1328	ssmin = ha/a;
1329	ssmax = fa*a;
1330	if( (double)(mm)==(double)(0) )
1331	{
1332
1333	//
1334	// Note that M is very tiny
1335	//
1336	if( (double)(l)==(double)(0) )
1337	{
1338	t = extsignbdsqr(2, ft)*extsignbdsqr(1, gt);
1339	}
1340	else
1341	{
1342	t = gt/extsignbdsqr(d, ft)+m/t;
1343	}
1344	}
1345	else
1346	{
1347	t = (m/(s+t)+m/(r+l))*(1+a);
1348	}
1349	l = Math.Sqrt(t*t+4);
1350	crt = 2/l;
1351	srt = t/l;
1352	clt = (crt+srt*m)/a;
1353	v = ht/ft;
1354	slt = v*srt/a;
1355	}
1356	}
1357	if( swp )
1358	{
1359	csl = srt;
1360	snl = crt;
1361	csr = slt;
1362	snr = clt;
1363	}
1364	else
1365	{
1366	csl = clt;
1367	snl = slt;
1368	csr = crt;
1369	snr = srt;
1370	}
1371
1372	//
1373	// Correct signs of SSMAX and SSMIN
1374	//
1375	if( pmax==1 )
1376	{
1377	tsign = extsignbdsqr(1, csr)extsignbdsqr(1, csl)extsignbdsqr(1, f);
1378	}
1379	if( pmax==2 )
1380	{
1381	tsign = extsignbdsqr(1, snr)extsignbdsqr(1, csl)extsignbdsqr(1, g);
1382	}
1383	if( pmax==3 )
1384	{
1385	tsign = extsignbdsqr(1, snr)extsignbdsqr(1, snl)extsignbdsqr(1, h);
1386	}
1387	ssmax = extsignbdsqr(ssmax, tsign);
1388	ssmin = extsignbdsqr(ssmin, tsignextsignbdsqr(1, f)extsignbdsqr(1, h));
1389	}
1390	}
1391	}

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

Download in other formats:

Update cookies preferences