Context Navigation

source: trunk/sources/ALGLIB/bidiagonal.cs @ 2474

Visit:

Last change on this file since 2474 was 2430, checked in by gkronber, 15 years ago
Moved ALGLIB code into a separate plugin. #783
File size: 48.5 KB

Rev	Line
[2154]	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
[2430]	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.
[2154]	15
[2430]	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.
[2154]	20
[2430]	21	A copy of the GNU General Public License is available at
	22	http://www.fsf.org/licensing/licenses
[2154]	23
[2430]	24	>>> END OF LICENSE >>>
[2154]	25	*************************************************************************/
	26
	27	using System;
	28
[2430]	29	namespace alglib
[2154]	30	{
[2430]	31	public class bidiagonal
	32	{
	33	/*************************************************************************
	34	Reduction of a rectangular matrix to bidiagonal form
[2154]	35
[2430]	36	The algorithm reduces the rectangular matrix A to bidiagonal form by
	37	orthogonal transformations P and Q: A = QBP.
[2154]	38
[2430]	39	Input parameters:
	40	A - source matrix. array[0..M-1, 0..N-1]
	41	M - number of rows in matrix A.
	42	N - number of columns in matrix A.
[2154]	43
[2430]	44	Output parameters:
	45	A - matrices Q, B, P in compact form (see below).
	46	TauQ - scalar factors which are used to form matrix Q.
	47	TauP - scalar factors which are used to form matrix P.
[2154]	48
[2430]	49	The main diagonal and one of the secondary diagonals of matrix A are
	50	replaced with bidiagonal matrix B. Other elements contain elementary
	51	reflections which form MxM matrix Q and NxN matrix P, respectively.
[2154]	52
[2430]	53	If M>=N, B is the upper bidiagonal MxN matrix and is stored in the
	54	corresponding elements of matrix A. Matrix Q is represented as a
	55	product of elementary reflections Q = H(0)H(1)...*H(n-1), where
	56	H(i) = 1-tauvv'. Here tau is a scalar which is stored in TauQ[i], and
	57	vector v has the following structure: v(0:i-1)=0, v(i)=1, v(i+1:m-1) is
	58	stored in elements A(i+1:m-1,i). Matrix P is as follows: P =
	59	G(0)G(1)...G(n-2), where G(i) = 1 - tauu*u'. Tau is stored in TauP[i],
	60	u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).
[2154]	61
[2430]	62	If M<N, B is the lower bidiagonal MxN matrix and is stored in the
	63	corresponding elements of matrix A. Q = H(0)H(1)...*H(m-2), where
	64	H(i) = 1 - tauvv', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
	65	is stored in elements A(i+2:m-1,i). P = G(0)G(1)...*G(m-1),
	66	G(i) = 1-tauuu', tau is stored in TauP, u(0:i-1)=0, u(i)=1, u(i+1:n-1)
	67	is stored in A(i,i+1:n-1).
[2154]	68
[2430]	69	EXAMPLE:
[2154]	70
[2430]	71	m=6, n=5 (m > n): m=5, n=6 (m < n):
[2154]	72
[2430]	73	( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
	74	( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
	75	( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
	76	( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
	77	( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
	78	( v1 v2 v3 v4 v5 )
[2154]	79
[2430]	80	Here vi and ui are vectors which form H(i) and G(i), and d and e -
	81	are the diagonal and off-diagonal elements of matrix B.
	82	*************************************************************************/
	83	public static void rmatrixbd(ref double[,] a,
	84	int m,
	85	int n,
	86	ref double[] tauq,
	87	ref double[] taup)
	88	{
	89	double[] work = new double[0];
	90	double[] t = new double[0];
	91	int minmn = 0;
	92	int maxmn = 0;
	93	int i = 0;
	94	int j = 0;
	95	double ltau = 0;
	96	int i_ = 0;
	97	int i1_ = 0;
[2154]	98
	99
	100	//
[2430]	101	// Prepare
[2154]	102	//
[2430]	103	if( n<=0 \| m<=0 )
[2154]	104	{
[2430]	105	return;
	106	}
	107	minmn = Math.Min(m, n);
	108	maxmn = Math.Max(m, n);
	109	work = new double[maxmn+1];
	110	t = new double[maxmn+1];
	111	if( m>=n )
	112	{
	113	tauq = new double[n-1+1];
	114	taup = new double[n-1+1];
	115	}
	116	else
	117	{
	118	tauq = new double[m-1+1];
	119	taup = new double[m-1+1];
	120	}
	121	if( m>=n )
	122	{
[2154]	123
	124	//
[2430]	125	// Reduce to upper bidiagonal form
[2154]	126	//
[2430]	127	for(i=0; i<=n-1; i++)
[2154]	128	{
	129
	130	//
[2430]	131	// Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
[2154]	132	//
[2430]	133	i1_ = (i) - (1);
	134	for(i_=1; i_<=m-i;i_++)
[2154]	135	{
[2430]	136	t[i_] = a[i_+i1_,i];
[2154]	137	}
[2430]	138	reflections.generatereflection(ref t, m-i, ref ltau);
	139	tauq[i] = ltau;
	140	i1_ = (1) - (i);
	141	for(i_=i; i_<=m-1;i_++)
[2154]	142	{
[2430]	143	a[i_,i] = t[i_+i1_];
[2154]	144	}
	145	t[1] = 1;
	146
	147	//
[2430]	148	// Apply H(i) to A(i:m-1,i+1:n-1) from the left
[2154]	149	//
[2430]	150	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m-1, i+1, n-1, ref work);
	151	if( i<n-1 )
	152	{
	153
	154	//
	155	// Generate elementary reflector G(i) to annihilate
	156	// A(i,i+2:n-1)
	157	//
	158	i1_ = (i+1) - (1);
	159	for(i_=1; i_<=n-i-1;i_++)
	160	{
	161	t[i_] = a[i,i_+i1_];
	162	}
	163	reflections.generatereflection(ref t, n-1-i, ref ltau);
	164	taup[i] = ltau;
	165	i1_ = (1) - (i+1);
	166	for(i_=i+1; i_<=n-1;i_++)
	167	{
	168	a[i,i_] = t[i_+i1_];
	169	}
	170	t[1] = 1;
	171
	172	//
	173	// Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
	174	//
	175	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
	176	}
	177	else
	178	{
	179	taup[i] = 0;
	180	}
[2154]	181	}
	182	}
[2430]	183	else
[2154]	184	{
	185
	186	//
[2430]	187	// Reduce to lower bidiagonal form
[2154]	188	//
[2430]	189	for(i=0; i<=m-1; i++)
[2154]	190	{
	191
	192	//
[2430]	193	// Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
[2154]	194	//
[2430]	195	i1_ = (i) - (1);
	196	for(i_=1; i_<=n-i;i_++)
[2154]	197	{
[2430]	198	t[i_] = a[i,i_+i1_];
[2154]	199	}
[2430]	200	reflections.generatereflection(ref t, n-i, ref ltau);
	201	taup[i] = ltau;
	202	i1_ = (1) - (i);
	203	for(i_=i; i_<=n-1;i_++)
[2154]	204	{
[2430]	205	a[i,i_] = t[i_+i1_];
[2154]	206	}
	207	t[1] = 1;
	208
	209	//
[2430]	210	// Apply G(i) to A(i+1:m-1,i:n-1) from the right
[2154]	211	//
[2430]	212	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i, n-1, ref work);
	213	if( i<m-1 )
	214	{
	215
	216	//
	217	// Generate elementary reflector H(i) to annihilate
	218	// A(i+2:m-1,i)
	219	//
	220	i1_ = (i+1) - (1);
	221	for(i_=1; i_<=m-1-i;i_++)
	222	{
	223	t[i_] = a[i_+i1_,i];
	224	}
	225	reflections.generatereflection(ref t, m-1-i, ref ltau);
	226	tauq[i] = ltau;
	227	i1_ = (1) - (i+1);
	228	for(i_=i+1; i_<=m-1;i_++)
	229	{
	230	a[i_,i] = t[i_+i1_];
	231	}
	232	t[1] = 1;
	233
	234	//
	235	// Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
	236	//
	237	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
	238	}
	239	else
	240	{
	241	tauq[i] = 0;
	242	}
[2154]	243	}
	244	}
	245	}
	246
	247
[2430]	248	/*************************************************************************
	249	Unpacking matrix Q which reduces a matrix to bidiagonal form.
[2154]	250
[2430]	251	Input parameters:
	252	QP - matrices Q and P in compact form.
	253	Output of ToBidiagonal subroutine.
	254	M - number of rows in matrix A.
	255	N - number of columns in matrix A.
	256	TAUQ - scalar factors which are used to form Q.
	257	Output of ToBidiagonal subroutine.
	258	QColumns - required number of columns in matrix Q.
	259	M>=QColumns>=0.
[2154]	260
[2430]	261	Output parameters:
	262	Q - first QColumns columns of matrix Q.
	263	Array[0..M-1, 0..QColumns-1]
	264	If QColumns=0, the array is not modified.
[2154]	265
[2430]	266	-- ALGLIB --
	267	Copyright 2005 by Bochkanov Sergey
	268	*************************************************************************/
	269	public static void rmatrixbdunpackq(ref double[,] qp,
	270	int m,
	271	int n,
	272	ref double[] tauq,
	273	int qcolumns,
	274	ref double[,] q)
	275	{
	276	int i = 0;
	277	int j = 0;
[2154]	278
[2430]	279	System.Diagnostics.Debug.Assert(qcolumns<=m, "RMatrixBDUnpackQ: QColumns>M!");
	280	System.Diagnostics.Debug.Assert(qcolumns>=0, "RMatrixBDUnpackQ: QColumns<0!");
	281	if( m==0 \| n==0 \| qcolumns==0 )
[2154]	282	{
[2430]	283	return;
	284	}
	285
	286	//
	287	// prepare Q
	288	//
	289	q = new double[m-1+1, qcolumns-1+1];
	290	for(i=0; i<=m-1; i++)
	291	{
	292	for(j=0; j<=qcolumns-1; j++)
[2154]	293	{
[2430]	294	if( i==j )
	295	{
	296	q[i,j] = 1;
	297	}
	298	else
	299	{
	300	q[i,j] = 0;
	301	}
[2154]	302	}
	303	}
[2430]	304
	305	//
	306	// Calculate
	307	//
	308	rmatrixbdmultiplybyq(ref qp, m, n, ref tauq, ref q, m, qcolumns, false, false);
[2154]	309	}
	310
	311
[2430]	312	/*************************************************************************
	313	Multiplication by matrix Q which reduces matrix A to bidiagonal form.
[2154]	314
[2430]	315	The algorithm allows pre- or post-multiply by Q or Q'.
[2154]	316
[2430]	317	Input parameters:
	318	QP - matrices Q and P in compact form.
	319	Output of ToBidiagonal subroutine.
	320	M - number of rows in matrix A.
	321	N - number of columns in matrix A.
	322	TAUQ - scalar factors which are used to form Q.
	323	Output of ToBidiagonal subroutine.
	324	Z - multiplied matrix.
	325	array[0..ZRows-1,0..ZColumns-1]
	326	ZRows - number of rows in matrix Z. If FromTheRight=False,
	327	ZRows=M, otherwise ZRows can be arbitrary.
	328	ZColumns - number of columns in matrix Z. If FromTheRight=True,
	329	ZColumns=M, otherwise ZColumns can be arbitrary.
	330	FromTheRight - pre- or post-multiply.
	331	DoTranspose - multiply by Q or Q'.
[2154]	332
[2430]	333	Output parameters:
	334	Z - product of Z and Q.
	335	Array[0..ZRows-1,0..ZColumns-1]
	336	If ZRows=0 or ZColumns=0, the array is not modified.
[2154]	337
[2430]	338	-- ALGLIB --
	339	Copyright 2005 by Bochkanov Sergey
	340	*************************************************************************/
	341	public static void rmatrixbdmultiplybyq(ref double[,] qp,
	342	int m,
	343	int n,
	344	ref double[] tauq,
	345	ref double[,] z,
	346	int zrows,
	347	int zcolumns,
	348	bool fromtheright,
	349	bool dotranspose)
	350	{
	351	int i = 0;
	352	int i1 = 0;
	353	int i2 = 0;
	354	int istep = 0;
	355	double[] v = new double[0];
	356	double[] work = new double[0];
	357	int mx = 0;
	358	int i_ = 0;
	359	int i1_ = 0;
[2154]	360
[2430]	361	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
[2154]	362	{
[2430]	363	return;
[2154]	364	}
[2430]	365	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m \| !fromtheright & zrows==m, "RMatrixBDMultiplyByQ: incorrect Z size!");
[2154]	366
	367	//
[2430]	368	// init
[2154]	369	//
[2430]	370	mx = Math.Max(m, n);
	371	mx = Math.Max(mx, zrows);
	372	mx = Math.Max(mx, zcolumns);
	373	v = new double[mx+1];
	374	work = new double[mx+1];
	375	if( m>=n )
[2154]	376	{
[2430]	377
	378	//
	379	// setup
	380	//
[2154]	381	if( fromtheright )
	382	{
[2430]	383	i1 = 0;
	384	i2 = n-1;
	385	istep = +1;
[2154]	386	}
	387	else
	388	{
[2430]	389	i1 = n-1;
	390	i2 = 0;
	391	istep = -1;
[2154]	392	}
[2430]	393	if( dotranspose )
	394	{
	395	i = i1;
	396	i1 = i2;
	397	i2 = i;
	398	istep = -istep;
	399	}
	400
	401	//
	402	// Process
	403	//
[2154]	404	i = i1;
	405	do
	406	{
[2430]	407	i1_ = (i) - (1);
	408	for(i_=1; i_<=m-i;i_++)
[2154]	409	{
	410	v[i_] = qp[i_+i1_,i];
	411	}
	412	v[1] = 1;
	413	if( fromtheright )
	414	{
[2430]	415	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i, m-1, ref work);
[2154]	416	}
	417	else
	418	{
[2430]	419	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m-1, 0, zcolumns-1, ref work);
[2154]	420	}
	421	i = i+istep;
	422	}
	423	while( i!=i2+istep );
	424	}
[2430]	425	else
[2154]	426	{
[2430]	427
	428	//
	429	// setup
	430	//
	431	if( fromtheright )
[2154]	432	{
[2430]	433	i1 = 0;
	434	i2 = m-2;
	435	istep = +1;
[2154]	436	}
	437	else
	438	{
[2430]	439	i1 = m-2;
	440	i2 = 0;
	441	istep = -1;
[2154]	442	}
[2430]	443	if( dotranspose )
	444	{
	445	i = i1;
	446	i1 = i2;
	447	i2 = i;
	448	istep = -istep;
	449	}
	450
	451	//
	452	// Process
	453	//
	454	if( m-1>0 )
	455	{
	456	i = i1;
	457	do
	458	{
	459	i1_ = (i+1) - (1);
	460	for(i_=1; i_<=m-i-1;i_++)
	461	{
	462	v[i_] = qp[i_+i1_,i];
	463	}
	464	v[1] = 1;
	465	if( fromtheright )
	466	{
	467	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i+1, m-1, ref work);
	468	}
	469	else
	470	{
	471	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m-1, 0, zcolumns-1, ref work);
	472	}
	473	i = i+istep;
	474	}
	475	while( i!=i2+istep );
	476	}
[2154]	477	}
	478	}
	479
	480
[2430]	481	/*************************************************************************
	482	Unpacking matrix P which reduces matrix A to bidiagonal form.
	483	The subroutine returns transposed matrix P.
[2154]	484
[2430]	485	Input parameters:
	486	QP - matrices Q and P in compact form.
	487	Output of ToBidiagonal subroutine.
	488	M - number of rows in matrix A.
	489	N - number of columns in matrix A.
	490	TAUP - scalar factors which are used to form P.
	491	Output of ToBidiagonal subroutine.
	492	PTRows - required number of rows of matrix P^T. N >= PTRows >= 0.
[2154]	493
[2430]	494	Output parameters:
	495	PT - first PTRows columns of matrix P^T
	496	Array[0..PTRows-1, 0..N-1]
	497	If PTRows=0, the array is not modified.
[2154]	498
[2430]	499	-- ALGLIB --
	500	Copyright 2005-2007 by Bochkanov Sergey
	501	*************************************************************************/
	502	public static void rmatrixbdunpackpt(ref double[,] qp,
	503	int m,
	504	int n,
	505	ref double[] taup,
	506	int ptrows,
	507	ref double[,] pt)
	508	{
	509	int i = 0;
	510	int j = 0;
[2154]	511
[2430]	512	System.Diagnostics.Debug.Assert(ptrows<=n, "RMatrixBDUnpackPT: PTRows>N!");
	513	System.Diagnostics.Debug.Assert(ptrows>=0, "RMatrixBDUnpackPT: PTRows<0!");
	514	if( m==0 \| n==0 \| ptrows==0 )
[2154]	515	{
[2430]	516	return;
[2154]	517	}
	518
	519	//
[2430]	520	// prepare PT
[2154]	521	//
[2430]	522	pt = new double[ptrows-1+1, n-1+1];
	523	for(i=0; i<=ptrows-1; i++)
[2154]	524	{
[2430]	525	for(j=0; j<=n-1; j++)
[2154]	526	{
[2430]	527	if( i==j )
[2154]	528	{
[2430]	529	pt[i,j] = 1;
[2154]	530	}
	531	else
	532	{
[2430]	533	pt[i,j] = 0;
[2154]	534	}
	535	}
	536	}
	537
	538	//
[2430]	539	// Calculate
[2154]	540	//
[2430]	541	rmatrixbdmultiplybyp(ref qp, m, n, ref taup, ref pt, ptrows, n, true, true);
	542	}
	543
	544
	545	/*************************************************************************
	546	Multiplication by matrix P which reduces matrix A to bidiagonal form.
	547
	548	The algorithm allows pre- or post-multiply by P or P'.
	549
	550	Input parameters:
	551	QP - matrices Q and P in compact form.
	552	Output of RMatrixBD subroutine.
	553	M - number of rows in matrix A.
	554	N - number of columns in matrix A.
	555	TAUP - scalar factors which are used to form P.
	556	Output of RMatrixBD subroutine.
	557	Z - multiplied matrix.
	558	Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
	559	ZRows - number of rows in matrix Z. If FromTheRight=False,
	560	ZRows=N, otherwise ZRows can be arbitrary.
	561	ZColumns - number of columns in matrix Z. If FromTheRight=True,
	562	ZColumns=N, otherwise ZColumns can be arbitrary.
	563	FromTheRight - pre- or post-multiply.
	564	DoTranspose - multiply by P or P'.
	565
	566	Output parameters:
	567	Z - product of Z and P.
	568	Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
	569	If ZRows=0 or ZColumns=0, the array is not modified.
	570
	571	-- ALGLIB --
	572	Copyright 2005-2007 by Bochkanov Sergey
	573	*************************************************************************/
	574	public static void rmatrixbdmultiplybyp(ref double[,] qp,
	575	int m,
	576	int n,
	577	ref double[] taup,
	578	ref double[,] z,
	579	int zrows,
	580	int zcolumns,
	581	bool fromtheright,
	582	bool dotranspose)
	583	{
	584	int i = 0;
	585	double[] v = new double[0];
	586	double[] work = new double[0];
	587	int mx = 0;
	588	int i1 = 0;
	589	int i2 = 0;
	590	int istep = 0;
	591	int i_ = 0;
	592	int i1_ = 0;
	593
	594	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
[2154]	595	{
[2430]	596	return;
[2154]	597	}
[2430]	598	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n \| !fromtheright & zrows==n, "RMatrixBDMultiplyByP: incorrect Z size!");
[2154]	599
	600	//
[2430]	601	// init
[2154]	602	//
[2430]	603	mx = Math.Max(m, n);
	604	mx = Math.Max(mx, zrows);
	605	mx = Math.Max(mx, zcolumns);
	606	v = new double[mx+1];
	607	work = new double[mx+1];
	608	v = new double[mx+1];
	609	work = new double[mx+1];
	610	if( m>=n )
[2154]	611	{
[2430]	612
	613	//
	614	// setup
	615	//
	616	if( fromtheright )
[2154]	617	{
[2430]	618	i1 = n-2;
	619	i2 = 0;
	620	istep = -1;
[2154]	621	}
[2430]	622	else
	623	{
	624	i1 = 0;
	625	i2 = n-2;
	626	istep = +1;
	627	}
	628	if( !dotranspose )
	629	{
	630	i = i1;
	631	i1 = i2;
	632	i2 = i;
	633	istep = -istep;
	634	}
	635
	636	//
	637	// Process
	638	//
	639	if( n-1>0 )
	640	{
	641	i = i1;
	642	do
	643	{
	644	i1_ = (i+1) - (1);
	645	for(i_=1; i_<=n-1-i;i_++)
	646	{
	647	v[i_] = qp[i,i_+i1_];
	648	}
	649	v[1] = 1;
	650	if( fromtheright )
	651	{
	652	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i+1, n-1, ref work);
	653	}
	654	else
	655	{
	656	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n-1, 0, zcolumns-1, ref work);
	657	}
	658	i = i+istep;
	659	}
	660	while( i!=i2+istep );
	661	}
	662	}
	663	else
	664	{
	665
	666	//
	667	// setup
	668	//
[2154]	669	if( fromtheright )
	670	{
[2430]	671	i1 = m-1;
	672	i2 = 0;
	673	istep = -1;
[2154]	674	}
	675	else
	676	{
[2430]	677	i1 = 0;
	678	i2 = m-1;
	679	istep = +1;
[2154]	680	}
[2430]	681	if( !dotranspose )
	682	{
	683	i = i1;
	684	i1 = i2;
	685	i2 = i;
	686	istep = -istep;
	687	}
	688
	689	//
	690	// Process
	691	//
	692	i = i1;
	693	do
	694	{
	695	i1_ = (i) - (1);
	696	for(i_=1; i_<=n-i;i_++)
	697	{
	698	v[i_] = qp[i,i_+i1_];
	699	}
	700	v[1] = 1;
	701	if( fromtheright )
	702	{
	703	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i, n-1, ref work);
	704	}
	705	else
	706	{
	707	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n-1, 0, zcolumns-1, ref work);
	708	}
	709	i = i+istep;
	710	}
	711	while( i!=i2+istep );
[2154]	712	}
	713	}
	714
	715
[2430]	716	/*************************************************************************
	717	Unpacking of the main and secondary diagonals of bidiagonal decomposition
	718	of matrix A.
[2154]	719
[2430]	720	Input parameters:
	721	B - output of RMatrixBD subroutine.
	722	M - number of rows in matrix B.
	723	N - number of columns in matrix B.
[2154]	724
[2430]	725	Output parameters:
	726	IsUpper - True, if the matrix is upper bidiagonal.
	727	otherwise IsUpper is False.
	728	D - the main diagonal.
	729	Array whose index ranges within [0..Min(M,N)-1].
	730	E - the secondary diagonal (upper or lower, depending on
	731	the value of IsUpper).
	732	Array index ranges within [0..Min(M,N)-1], the last
	733	element is not used.
[2154]	734
[2430]	735	-- ALGLIB --
	736	Copyright 2005-2007 by Bochkanov Sergey
	737	*************************************************************************/
	738	public static void rmatrixbdunpackdiagonals(ref double[,] b,
	739	int m,
	740	int n,
	741	ref bool isupper,
	742	ref double[] d,
	743	ref double[] e)
	744	{
	745	int i = 0;
[2154]	746
[2430]	747	isupper = m>=n;
	748	if( m<=0 \| n<=0 )
[2154]	749	{
[2430]	750	return;
[2154]	751	}
[2430]	752	if( isupper )
[2154]	753	{
[2430]	754	d = new double[n-1+1];
	755	e = new double[n-1+1];
	756	for(i=0; i<=n-2; i++)
	757	{
	758	d[i] = b[i,i];
	759	e[i] = b[i,i+1];
	760	}
	761	d[n-1] = b[n-1,n-1];
[2154]	762	}
[2430]	763	else
	764	{
	765	d = new double[m-1+1];
	766	e = new double[m-1+1];
	767	for(i=0; i<=m-2; i++)
	768	{
	769	d[i] = b[i,i];
	770	e[i] = b[i+1,i];
	771	}
	772	d[m-1] = b[m-1,m-1];
	773	}
[2154]	774	}
	775
	776
[2430]	777	/*************************************************************************
	778	Obsolete 1-based subroutine.
	779	See RMatrixBD for 0-based replacement.
	780	*************************************************************************/
	781	public static void tobidiagonal(ref double[,] a,
	782	int m,
	783	int n,
	784	ref double[] tauq,
	785	ref double[] taup)
	786	{
	787	double[] work = new double[0];
	788	double[] t = new double[0];
	789	int minmn = 0;
	790	int maxmn = 0;
	791	int i = 0;
	792	double ltau = 0;
	793	int mmip1 = 0;
	794	int nmi = 0;
	795	int ip1 = 0;
	796	int nmip1 = 0;
	797	int mmi = 0;
	798	int i_ = 0;
	799	int i1_ = 0;
[2154]	800
[2430]	801	minmn = Math.Min(m, n);
	802	maxmn = Math.Max(m, n);
	803	work = new double[maxmn+1];
	804	t = new double[maxmn+1];
	805	taup = new double[minmn+1];
	806	tauq = new double[minmn+1];
	807	if( m>=n )
[2154]	808	{
	809
	810	//
[2430]	811	// Reduce to upper bidiagonal form
[2154]	812	//
[2430]	813	for(i=1; i<=n; i++)
[2154]	814	{
	815
	816	//
[2430]	817	// Generate elementary reflector H(i) to annihilate A(i+1:m,i)
[2154]	818	//
[2430]	819	mmip1 = m-i+1;
	820	i1_ = (i) - (1);
	821	for(i_=1; i_<=mmip1;i_++)
[2154]	822	{
[2430]	823	t[i_] = a[i_+i1_,i];
[2154]	824	}
[2430]	825	reflections.generatereflection(ref t, mmip1, ref ltau);
	826	tauq[i] = ltau;
	827	i1_ = (1) - (i);
	828	for(i_=i; i_<=m;i_++)
[2154]	829	{
[2430]	830	a[i_,i] = t[i_+i1_];
[2154]	831	}
	832	t[1] = 1;
	833
	834	//
[2430]	835	// Apply H(i) to A(i:m,i+1:n) from the left
[2154]	836	//
[2430]	837	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m, i+1, n, ref work);
	838	if( i<n )
	839	{
	840
	841	//
	842	// Generate elementary reflector G(i) to annihilate
	843	// A(i,i+2:n)
	844	//
	845	nmi = n-i;
	846	ip1 = i+1;
	847	i1_ = (ip1) - (1);
	848	for(i_=1; i_<=nmi;i_++)
	849	{
	850	t[i_] = a[i,i_+i1_];
	851	}
	852	reflections.generatereflection(ref t, nmi, ref ltau);
	853	taup[i] = ltau;
	854	i1_ = (1) - (ip1);
	855	for(i_=ip1; i_<=n;i_++)
	856	{
	857	a[i,i_] = t[i_+i1_];
	858	}
	859	t[1] = 1;
	860
	861	//
	862	// Apply G(i) to A(i+1:m,i+1:n) from the right
	863	//
	864	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
	865	}
	866	else
	867	{
	868	taup[i] = 0;
	869	}
[2154]	870	}
	871	}
[2430]	872	else
[2154]	873	{
	874
	875	//
[2430]	876	// Reduce to lower bidiagonal form
[2154]	877	//
[2430]	878	for(i=1; i<=m; i++)
[2154]	879	{
	880
	881	//
[2430]	882	// Generate elementary reflector G(i) to annihilate A(i,i+1:n)
[2154]	883	//
[2430]	884	nmip1 = n-i+1;
	885	i1_ = (i) - (1);
	886	for(i_=1; i_<=nmip1;i_++)
[2154]	887	{
[2430]	888	t[i_] = a[i,i_+i1_];
[2154]	889	}
[2430]	890	reflections.generatereflection(ref t, nmip1, ref ltau);
	891	taup[i] = ltau;
	892	i1_ = (1) - (i);
	893	for(i_=i; i_<=n;i_++)
[2154]	894	{
[2430]	895	a[i,i_] = t[i_+i1_];
[2154]	896	}
	897	t[1] = 1;
	898
	899	//
[2430]	900	// Apply G(i) to A(i+1:m,i:n) from the right
[2154]	901	//
[2430]	902	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m, i, n, ref work);
	903	if( i<m )
	904	{
	905
	906	//
	907	// Generate elementary reflector H(i) to annihilate
	908	// A(i+2:m,i)
	909	//
	910	mmi = m-i;
	911	ip1 = i+1;
	912	i1_ = (ip1) - (1);
	913	for(i_=1; i_<=mmi;i_++)
	914	{
	915	t[i_] = a[i_+i1_,i];
	916	}
	917	reflections.generatereflection(ref t, mmi, ref ltau);
	918	tauq[i] = ltau;
	919	i1_ = (1) - (ip1);
	920	for(i_=ip1; i_<=m;i_++)
	921	{
	922	a[i_,i] = t[i_+i1_];
	923	}
	924	t[1] = 1;
	925
	926	//
	927	// Apply H(i) to A(i+1:m,i+1:n) from the left
	928	//
	929	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m, i+1, n, ref work);
	930	}
	931	else
	932	{
	933	tauq[i] = 0;
	934	}
[2154]	935	}
	936	}
	937	}
	938
	939
[2430]	940	/*************************************************************************
	941	Obsolete 1-based subroutine.
	942	See RMatrixBDUnpackQ for 0-based replacement.
	943	*************************************************************************/
	944	public static void unpackqfrombidiagonal(ref double[,] qp,
	945	int m,
	946	int n,
	947	ref double[] tauq,
	948	int qcolumns,
	949	ref double[,] q)
	950	{
	951	int i = 0;
	952	int j = 0;
	953	int ip1 = 0;
	954	double[] v = new double[0];
	955	double[] work = new double[0];
	956	int vm = 0;
	957	int i_ = 0;
	958	int i1_ = 0;
[2154]	959
[2430]	960	System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromBidiagonal: QColumns>M!");
	961	if( m==0 \| n==0 \| qcolumns==0 )
[2154]	962	{
[2430]	963	return;
	964	}
	965
	966	//
	967	// init
	968	//
	969	q = new double[m+1, qcolumns+1];
	970	v = new double[m+1];
	971	work = new double[qcolumns+1];
	972
	973	//
	974	// prepare Q
	975	//
	976	for(i=1; i<=m; i++)
	977	{
	978	for(j=1; j<=qcolumns; j++)
[2154]	979	{
[2430]	980	if( i==j )
	981	{
	982	q[i,j] = 1;
	983	}
	984	else
	985	{
	986	q[i,j] = 0;
	987	}
[2154]	988	}
	989	}
[2430]	990	if( m>=n )
[2154]	991	{
[2430]	992	for(i=Math.Min(n, qcolumns); i>=1; i--)
[2154]	993	{
[2430]	994	vm = m-i+1;
	995	i1_ = (i) - (1);
	996	for(i_=1; i_<=vm;i_++)
	997	{
	998	v[i_] = qp[i_+i1_,i];
	999	}
	1000	v[1] = 1;
	1001	reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i, m, 1, qcolumns, ref work);
[2154]	1002	}
	1003	}
[2430]	1004	else
[2154]	1005	{
[2430]	1006	for(i=Math.Min(m-1, qcolumns-1); i>=1; i--)
[2154]	1007	{
[2430]	1008	vm = m-i;
	1009	ip1 = i+1;
	1010	i1_ = (ip1) - (1);
	1011	for(i_=1; i_<=vm;i_++)
	1012	{
	1013	v[i_] = qp[i_+i1_,i];
	1014	}
	1015	v[1] = 1;
	1016	reflections.applyreflectionfromtheleft(ref q, tauq[i], ref v, i+1, m, 1, qcolumns, ref work);
[2154]	1017	}
	1018	}
	1019	}
	1020
	1021
[2430]	1022	/*************************************************************************
	1023	Obsolete 1-based subroutine.
	1024	See RMatrixBDMultiplyByQ for 0-based replacement.
	1025	*************************************************************************/
	1026	public static void multiplybyqfrombidiagonal(ref double[,] qp,
	1027	int m,
	1028	int n,
	1029	ref double[] tauq,
	1030	ref double[,] z,
	1031	int zrows,
	1032	int zcolumns,
	1033	bool fromtheright,
	1034	bool dotranspose)
	1035	{
	1036	int i = 0;
	1037	int ip1 = 0;
	1038	int i1 = 0;
	1039	int i2 = 0;
	1040	int istep = 0;
	1041	double[] v = new double[0];
	1042	double[] work = new double[0];
	1043	int vm = 0;
	1044	int mx = 0;
	1045	int i_ = 0;
	1046	int i1_ = 0;
[2154]	1047
[2430]	1048	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
[2154]	1049	{
[2430]	1050	return;
[2154]	1051	}
[2430]	1052	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m \| !fromtheright & zrows==m, "MultiplyByQFromBidiagonal: incorrect Z size!");
[2154]	1053
	1054	//
[2430]	1055	// init
[2154]	1056	//
[2430]	1057	mx = Math.Max(m, n);
	1058	mx = Math.Max(mx, zrows);
	1059	mx = Math.Max(mx, zcolumns);
	1060	v = new double[mx+1];
	1061	work = new double[mx+1];
	1062	if( m>=n )
[2154]	1063	{
[2430]	1064
	1065	//
	1066	// setup
	1067	//
[2154]	1068	if( fromtheright )
	1069	{
[2430]	1070	i1 = 1;
	1071	i2 = n;
	1072	istep = +1;
[2154]	1073	}
	1074	else
	1075	{
[2430]	1076	i1 = n;
	1077	i2 = 1;
	1078	istep = -1;
[2154]	1079	}
[2430]	1080	if( dotranspose )
	1081	{
	1082	i = i1;
	1083	i1 = i2;
	1084	i2 = i;
	1085	istep = -istep;
	1086	}
	1087
	1088	//
	1089	// Process
	1090	//
[2154]	1091	i = i1;
	1092	do
	1093	{
[2430]	1094	vm = m-i+1;
	1095	i1_ = (i) - (1);
[2154]	1096	for(i_=1; i_<=vm;i_++)
	1097	{
	1098	v[i_] = qp[i_+i1_,i];
	1099	}
	1100	v[1] = 1;
	1101	if( fromtheright )
	1102	{
[2430]	1103	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i, m, ref work);
[2154]	1104	}
	1105	else
	1106	{
[2430]	1107	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m, 1, zcolumns, ref work);
[2154]	1108	}
	1109	i = i+istep;
	1110	}
	1111	while( i!=i2+istep );
	1112	}
[2430]	1113	else
[2154]	1114	{
[2430]	1115
	1116	//
	1117	// setup
	1118	//
	1119	if( fromtheright )
[2154]	1120	{
[2430]	1121	i1 = 1;
	1122	i2 = m-1;
	1123	istep = +1;
[2154]	1124	}
	1125	else
	1126	{
[2430]	1127	i1 = m-1;
	1128	i2 = 1;
	1129	istep = -1;
[2154]	1130	}
[2430]	1131	if( dotranspose )
[2154]	1132	{
[2430]	1133	i = i1;
	1134	i1 = i2;
	1135	i2 = i;
	1136	istep = -istep;
[2154]	1137	}
[2430]	1138
	1139	//
	1140	// Process
	1141	//
	1142	if( m-1>0 )
[2154]	1143	{
[2430]	1144	i = i1;
	1145	do
	1146	{
	1147	vm = m-i;
	1148	ip1 = i+1;
	1149	i1_ = (ip1) - (1);
	1150	for(i_=1; i_<=vm;i_++)
	1151	{
	1152	v[i_] = qp[i_+i1_,i];
	1153	}
	1154	v[1] = 1;
	1155	if( fromtheright )
	1156	{
	1157	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 1, zrows, i+1, m, ref work);
	1158	}
	1159	else
	1160	{
	1161	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m, 1, zcolumns, ref work);
	1162	}
	1163	i = i+istep;
	1164	}
	1165	while( i!=i2+istep );
[2154]	1166	}
	1167	}
	1168	}
	1169
	1170
[2430]	1171	/*************************************************************************
	1172	Obsolete 1-based subroutine.
	1173	See RMatrixBDUnpackPT for 0-based replacement.
	1174	*************************************************************************/
	1175	public static void unpackptfrombidiagonal(ref double[,] qp,
	1176	int m,
	1177	int n,
	1178	ref double[] taup,
	1179	int ptrows,
	1180	ref double[,] pt)
	1181	{
	1182	int i = 0;
	1183	int j = 0;
	1184	int ip1 = 0;
	1185	double[] v = new double[0];
	1186	double[] work = new double[0];
	1187	int vm = 0;
	1188	int i_ = 0;
	1189	int i1_ = 0;
[2154]	1190
[2430]	1191	System.Diagnostics.Debug.Assert(ptrows<=n, "UnpackPTFromBidiagonal: PTRows>N!");
	1192	if( m==0 \| n==0 \| ptrows==0 )
	1193	{
	1194	return;
	1195	}
[2154]	1196
	1197	//
[2430]	1198	// init
[2154]	1199	//
[2430]	1200	pt = new double[ptrows+1, n+1];
	1201	v = new double[n+1];
	1202	work = new double[ptrows+1];
[2154]	1203
	1204	//
[2430]	1205	// prepare PT
[2154]	1206	//
[2430]	1207	for(i=1; i<=ptrows; i++)
[2154]	1208	{
[2430]	1209	for(j=1; j<=n; j++)
[2154]	1210	{
[2430]	1211	if( i==j )
	1212	{
	1213	pt[i,j] = 1;
	1214	}
	1215	else
	1216	{
	1217	pt[i,j] = 0;
	1218	}
	1219	}
	1220	}
	1221	if( m>=n )
	1222	{
	1223	for(i=Math.Min(n-1, ptrows-1); i>=1; i--)
	1224	{
[2154]	1225	vm = n-i;
	1226	ip1 = i+1;
	1227	i1_ = (ip1) - (1);
	1228	for(i_=1; i_<=vm;i_++)
	1229	{
	1230	v[i_] = qp[i,i_+i1_];
	1231	}
	1232	v[1] = 1;
[2430]	1233	reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i+1, n, ref work);
	1234	}
	1235	}
	1236	else
	1237	{
	1238	for(i=Math.Min(m, ptrows); i>=1; i--)
	1239	{
	1240	vm = n-i+1;
	1241	i1_ = (i) - (1);
	1242	for(i_=1; i_<=vm;i_++)
[2154]	1243	{
[2430]	1244	v[i_] = qp[i,i_+i1_];
[2154]	1245	}
[2430]	1246	v[1] = 1;
	1247	reflections.applyreflectionfromtheright(ref pt, taup[i], ref v, 1, ptrows, i, n, ref work);
[2154]	1248	}
	1249	}
	1250	}
[2430]	1251
	1252
	1253	/*************************************************************************
	1254	Obsolete 1-based subroutine.
	1255	See RMatrixBDMultiplyByP for 0-based replacement.
	1256	*************************************************************************/
	1257	public static void multiplybypfrombidiagonal(ref double[,] qp,
	1258	int m,
	1259	int n,
	1260	ref double[] taup,
	1261	ref double[,] z,
	1262	int zrows,
	1263	int zcolumns,
	1264	bool fromtheright,
	1265	bool dotranspose)
[2154]	1266	{
[2430]	1267	int i = 0;
	1268	int ip1 = 0;
	1269	double[] v = new double[0];
	1270	double[] work = new double[0];
	1271	int vm = 0;
	1272	int mx = 0;
	1273	int i1 = 0;
	1274	int i2 = 0;
	1275	int istep = 0;
	1276	int i_ = 0;
	1277	int i1_ = 0;
	1278
	1279	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
[2154]	1280	{
[2430]	1281	return;
[2154]	1282	}
[2430]	1283	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n \| !fromtheright & zrows==n, "MultiplyByQFromBidiagonal: incorrect Z size!");
[2154]	1284
	1285	//
[2430]	1286	// init
[2154]	1287	//
[2430]	1288	mx = Math.Max(m, n);
	1289	mx = Math.Max(mx, zrows);
	1290	mx = Math.Max(mx, zcolumns);
	1291	v = new double[mx+1];
	1292	work = new double[mx+1];
	1293	v = new double[mx+1];
	1294	work = new double[mx+1];
	1295	if( m>=n )
[2154]	1296	{
[2430]	1297
	1298	//
	1299	// setup
	1300	//
	1301	if( fromtheright )
[2154]	1302	{
[2430]	1303	i1 = n-1;
	1304	i2 = 1;
	1305	istep = -1;
[2154]	1306	}
[2430]	1307	else
	1308	{
	1309	i1 = 1;
	1310	i2 = n-1;
	1311	istep = +1;
	1312	}
	1313	if( !dotranspose )
	1314	{
	1315	i = i1;
	1316	i1 = i2;
	1317	i2 = i;
	1318	istep = -istep;
	1319	}
	1320
	1321	//
	1322	// Process
	1323	//
	1324	if( n-1>0 )
	1325	{
	1326	i = i1;
	1327	do
	1328	{
	1329	vm = n-i;
	1330	ip1 = i+1;
	1331	i1_ = (ip1) - (1);
	1332	for(i_=1; i_<=vm;i_++)
	1333	{
	1334	v[i_] = qp[i,i_+i1_];
	1335	}
	1336	v[1] = 1;
	1337	if( fromtheright )
	1338	{
	1339	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i+1, n, ref work);
	1340	}
	1341	else
	1342	{
	1343	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n, 1, zcolumns, ref work);
	1344	}
	1345	i = i+istep;
	1346	}
	1347	while( i!=i2+istep );
	1348	}
	1349	}
	1350	else
	1351	{
	1352
	1353	//
	1354	// setup
	1355	//
[2154]	1356	if( fromtheright )
	1357	{
[2430]	1358	i1 = m;
	1359	i2 = 1;
	1360	istep = -1;
[2154]	1361	}
	1362	else
	1363	{
[2430]	1364	i1 = 1;
	1365	i2 = m;
	1366	istep = +1;
[2154]	1367	}
[2430]	1368	if( !dotranspose )
	1369	{
	1370	i = i1;
	1371	i1 = i2;
	1372	i2 = i;
	1373	istep = -istep;
	1374	}
	1375
	1376	//
	1377	// Process
	1378	//
	1379	i = i1;
	1380	do
	1381	{
	1382	vm = n-i+1;
	1383	i1_ = (i) - (1);
	1384	for(i_=1; i_<=vm;i_++)
	1385	{
	1386	v[i_] = qp[i,i_+i1_];
	1387	}
	1388	v[1] = 1;
	1389	if( fromtheright )
	1390	{
	1391	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 1, zrows, i, n, ref work);
	1392	}
	1393	else
	1394	{
	1395	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n, 1, zcolumns, ref work);
	1396	}
	1397	i = i+istep;
	1398	}
	1399	while( i!=i2+istep );
[2154]	1400	}
	1401	}
	1402
	1403
[2430]	1404	/*************************************************************************
	1405	Obsolete 1-based subroutine.
	1406	See RMatrixBDUnpackDiagonals for 0-based replacement.
	1407	*************************************************************************/
	1408	public static void unpackdiagonalsfrombidiagonal(ref double[,] b,
	1409	int m,
	1410	int n,
	1411	ref bool isupper,
	1412	ref double[] d,
	1413	ref double[] e)
	1414	{
	1415	int i = 0;
[2154]	1416
[2430]	1417	isupper = m>=n;
	1418	if( m==0 \| n==0 )
[2154]	1419	{
[2430]	1420	return;
[2154]	1421	}
[2430]	1422	if( isupper )
[2154]	1423	{
[2430]	1424	d = new double[n+1];
	1425	e = new double[n+1];
	1426	for(i=1; i<=n-1; i++)
	1427	{
	1428	d[i] = b[i,i];
	1429	e[i] = b[i,i+1];
	1430	}
	1431	d[n] = b[n,n];
[2154]	1432	}
[2430]	1433	else
	1434	{
	1435	d = new double[m+1];
	1436	e = new double[m+1];
	1437	for(i=1; i<=m-1; i++)
	1438	{
	1439	d[i] = b[i,i];
	1440	e[i] = b[i+1,i];
	1441	}
	1442	d[m] = b[m,m];
	1443	}
[2154]	1444	}
	1445	}
	1446	}

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

Download in other formats:

Update cookies preferences