Context Navigation

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/ortfac.cs @ 5192

Visit:

Last change on this file since 5192 was 3839, checked in by mkommend, 14 years ago
implemented first version of LR (ticket #1012)
File size: 132.5 KB

Line
1	/*************************************************************************
2	Copyright (c) 2005-2010 Sergey Bochkanov.
3
4	Additional copyrights:
5	1992-2007 The University of Tennessee (as indicated in subroutines
6	comments).
7
8	>>> SOURCE LICENSE >>>
9	This program is free software; you can redistribute it and/or modify
10	it under the terms of the GNU General Public License as published by
11	the Free Software Foundation (www.fsf.org); either version 2 of the
12	License, or (at your option) any later version.
13
14	This program is distributed in the hope that it will be useful,
15	but WITHOUT ANY WARRANTY; without even the implied warranty of
16	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17	GNU General Public License for more details.
18
19	A copy of the GNU General Public License is available at
20	http://www.fsf.org/licensing/licenses
21
22	>>> END OF LICENSE >>>
23	*************************************************************************/
24
25	using System;
26
27	namespace alglib
28	{
29	public class ortfac
30	{
31	/*************************************************************************
32	QR decomposition of a rectangular matrix of size MxN
33
34	Input parameters:
35	A - matrix A whose indexes range within [0..M-1, 0..N-1].
36	M - number of rows in matrix A.
37	N - number of columns in matrix A.
38
39	Output parameters:
40	A - matrices Q and R in compact form (see below).
41	Tau - array of scalar factors which are used to form
42	matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].
43
44	Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
45	MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.
46
47	The elements of matrix R are located on and above the main diagonal of
48	matrix A. The elements which are located in Tau array and below the main
49	diagonal of matrix A are used to form matrix Q as follows:
50
51	Matrix Q is represented as a product of elementary reflections
52
53	Q = H(0)H(2)...*H(k-1),
54
55	where k = min(m,n), and each H(i) is in the form
56
57	H(i) = 1 - tau * v * (v^T)
58
59	where tau is a scalar stored in Tau[I]; v - real vector,
60	so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).
61
62	-- ALGLIB routine --
63	17.02.2010
64	Bochkanov Sergey
65	*************************************************************************/
66	public static void rmatrixqr(ref double[,] a,
67	int m,
68	int n,
69	ref double[] tau)
70	{
71	double[] work = new double[0];
72	double[] t = new double[0];
73	double[] taubuf = new double[0];
74	int minmn = 0;
75	double[,] tmpa = new double[0,0];
76	double[,] tmpt = new double[0,0];
77	double[,] tmpr = new double[0,0];
78	int blockstart = 0;
79	int blocksize = 0;
80	int rowscount = 0;
81	int i = 0;
82	int j = 0;
83	int k = 0;
84	double v = 0;
85	int i_ = 0;
86	int i1_ = 0;
87
88	if( m<=0 \| n<=0 )
89	{
90	return;
91	}
92	minmn = Math.Min(m, n);
93	work = new double[Math.Max(m, n)+1];
94	t = new double[Math.Max(m, n)+1];
95	tau = new double[minmn];
96	taubuf = new double[minmn];
97	tmpa = new double[m, ablas.ablasblocksize(ref a)];
98	tmpt = new double[ablas.ablasblocksize(ref a), 2*ablas.ablasblocksize(ref a)];
99	tmpr = new double[2*ablas.ablasblocksize(ref a), n];
100
101	//
102	// Blocked code
103	//
104	blockstart = 0;
105	while( blockstart!=minmn )
106	{
107
108	//
109	// Determine block size
110	//
111	blocksize = minmn-blockstart;
112	if( blocksize>ablas.ablasblocksize(ref a) )
113	{
114	blocksize = ablas.ablasblocksize(ref a);
115	}
116	rowscount = m-blockstart;
117
118	//
119	// QR decomposition of submatrix.
120	// Matrix is copied to temporary storage to solve
121	// some TLB issues arising from non-contiguous memory
122	// access pattern.
123	//
124	ablas.rmatrixcopy(rowscount, blocksize, ref a, blockstart, blockstart, ref tmpa, 0, 0);
125	rmatrixqrbasecase(ref tmpa, rowscount, blocksize, ref work, ref t, ref taubuf);
126	ablas.rmatrixcopy(rowscount, blocksize, ref tmpa, 0, 0, ref a, blockstart, blockstart);
127	i1_ = (0) - (blockstart);
128	for(i_=blockstart; i_<=blockstart+blocksize-1;i_++)
129	{
130	tau[i_] = taubuf[i_+i1_];
131	}
132
133	//
134	// Update the rest, choose between:
135	// a) Level 2 algorithm (when the rest of the matrix is small enough)
136	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
137	// representation for products of Householder transformations',
138	// by R. Schreiber and C. Van Loan.
139	//
140	if( blockstart+blocksize<=n-1 )
141	{
142	if( n-blockstart-blocksize>=2ablas.ablasblocksize(ref a) \| rowscount>=4ablas.ablasblocksize(ref a) )
143	{
144
145	//
146	// Prepare block reflector
147	//
148	rmatrixblockreflector(ref tmpa, ref taubuf, true, rowscount, blocksize, ref tmpt, ref work);
149
150	//
151	// Multiply the rest of A by Q'.
152	//
153	// Q = E + YTY' = E + TmpATmpTTmpA'
154	// Q' = E + YT'Y' = E + TmpATmpT'TmpA'
155	//
156	ablas.rmatrixgemm(blocksize, n-blockstart-blocksize, rowscount, 1.0, ref tmpa, 0, 0, 1, ref a, blockstart, blockstart+blocksize, 0, 0.0, ref tmpr, 0, 0);
157	ablas.rmatrixgemm(blocksize, n-blockstart-blocksize, blocksize, 1.0, ref tmpt, 0, 0, 1, ref tmpr, 0, 0, 0, 0.0, ref tmpr, blocksize, 0);
158	ablas.rmatrixgemm(rowscount, n-blockstart-blocksize, blocksize, 1.0, ref tmpa, 0, 0, 0, ref tmpr, blocksize, 0, 0, 1.0, ref a, blockstart, blockstart+blocksize);
159	}
160	else
161	{
162
163	//
164	// Level 2 algorithm
165	//
166	for(i=0; i<=blocksize-1; i++)
167	{
168	i1_ = (i) - (1);
169	for(i_=1; i_<=rowscount-i;i_++)
170	{
171	t[i_] = tmpa[i_+i1_,i];
172	}
173	t[1] = 1;
174	reflections.applyreflectionfromtheleft(ref a, taubuf[i], ref t, blockstart+i, m-1, blockstart+blocksize, n-1, ref work);
175	}
176	}
177	}
178
179	//
180	// Advance
181	//
182	blockstart = blockstart+blocksize;
183	}
184	}
185
186
187	/*************************************************************************
188	LQ decomposition of a rectangular matrix of size MxN
189
190	Input parameters:
191	A - matrix A whose indexes range within [0..M-1, 0..N-1].
192	M - number of rows in matrix A.
193	N - number of columns in matrix A.
194
195	Output parameters:
196	A - matrices L and Q in compact form (see below)
197	Tau - array of scalar factors which are used to form
198	matrix Q. Array whose index ranges within [0..Min(M,N)-1].
199
200	Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size
201	MxM, L - lower triangular (or lower trapezoid) matrix of size M x N.
202
203	The elements of matrix L are located on and below the main diagonal of
204	matrix A. The elements which are located in Tau array and above the main
205	diagonal of matrix A are used to form matrix Q as follows:
206
207	Matrix Q is represented as a product of elementary reflections
208
209	Q = H(k-1)H(k-2)...H(1)H(0),
210
211	where k = min(m,n), and each H(i) is of the form
212
213	H(i) = 1 - tau * v * (v^T)
214
215	where tau is a scalar stored in Tau[I]; v - real vector, so that v(0:i-1)=0,
216	v(i) = 1, v(i+1:n-1) stored in A(i,i+1:n-1).
217
218	-- ALGLIB routine --
219	17.02.2010
220	Bochkanov Sergey
221	*************************************************************************/
222	public static void rmatrixlq(ref double[,] a,
223	int m,
224	int n,
225	ref double[] tau)
226	{
227	double[] work = new double[0];
228	double[] t = new double[0];
229	double[] taubuf = new double[0];
230	int minmn = 0;
231	double[,] tmpa = new double[0,0];
232	double[,] tmpt = new double[0,0];
233	double[,] tmpr = new double[0,0];
234	int blockstart = 0;
235	int blocksize = 0;
236	int columnscount = 0;
237	int i = 0;
238	int j = 0;
239	int k = 0;
240	double v = 0;
241	int i_ = 0;
242	int i1_ = 0;
243
244	if( m<=0 \| n<=0 )
245	{
246	return;
247	}
248	minmn = Math.Min(m, n);
249	work = new double[Math.Max(m, n)+1];
250	t = new double[Math.Max(m, n)+1];
251	tau = new double[minmn];
252	taubuf = new double[minmn];
253	tmpa = new double[ablas.ablasblocksize(ref a), n];
254	tmpt = new double[ablas.ablasblocksize(ref a), 2*ablas.ablasblocksize(ref a)];
255	tmpr = new double[m, 2*ablas.ablasblocksize(ref a)];
256
257	//
258	// Blocked code
259	//
260	blockstart = 0;
261	while( blockstart!=minmn )
262	{
263
264	//
265	// Determine block size
266	//
267	blocksize = minmn-blockstart;
268	if( blocksize>ablas.ablasblocksize(ref a) )
269	{
270	blocksize = ablas.ablasblocksize(ref a);
271	}
272	columnscount = n-blockstart;
273
274	//
275	// LQ decomposition of submatrix.
276	// Matrix is copied to temporary storage to solve
277	// some TLB issues arising from non-contiguous memory
278	// access pattern.
279	//
280	ablas.rmatrixcopy(blocksize, columnscount, ref a, blockstart, blockstart, ref tmpa, 0, 0);
281	rmatrixlqbasecase(ref tmpa, blocksize, columnscount, ref work, ref t, ref taubuf);
282	ablas.rmatrixcopy(blocksize, columnscount, ref tmpa, 0, 0, ref a, blockstart, blockstart);
283	i1_ = (0) - (blockstart);
284	for(i_=blockstart; i_<=blockstart+blocksize-1;i_++)
285	{
286	tau[i_] = taubuf[i_+i1_];
287	}
288
289	//
290	// Update the rest, choose between:
291	// a) Level 2 algorithm (when the rest of the matrix is small enough)
292	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
293	// representation for products of Householder transformations',
294	// by R. Schreiber and C. Van Loan.
295	//
296	if( blockstart+blocksize<=m-1 )
297	{
298	if( m-blockstart-blocksize>=2*ablas.ablasblocksize(ref a) )
299	{
300
301	//
302	// Prepare block reflector
303	//
304	rmatrixblockreflector(ref tmpa, ref taubuf, false, columnscount, blocksize, ref tmpt, ref work);
305
306	//
307	// Multiply the rest of A by Q.
308	//
309	// Q = E + YTY' = E + TmpA'TmpTTmpA
310	//
311	ablas.rmatrixgemm(m-blockstart-blocksize, blocksize, columnscount, 1.0, ref a, blockstart+blocksize, blockstart, 0, ref tmpa, 0, 0, 1, 0.0, ref tmpr, 0, 0);
312	ablas.rmatrixgemm(m-blockstart-blocksize, blocksize, blocksize, 1.0, ref tmpr, 0, 0, 0, ref tmpt, 0, 0, 0, 0.0, ref tmpr, 0, blocksize);
313	ablas.rmatrixgemm(m-blockstart-blocksize, columnscount, blocksize, 1.0, ref tmpr, 0, blocksize, 0, ref tmpa, 0, 0, 0, 1.0, ref a, blockstart+blocksize, blockstart);
314	}
315	else
316	{
317
318	//
319	// Level 2 algorithm
320	//
321	for(i=0; i<=blocksize-1; i++)
322	{
323	i1_ = (i) - (1);
324	for(i_=1; i_<=columnscount-i;i_++)
325	{
326	t[i_] = tmpa[i,i_+i1_];
327	}
328	t[1] = 1;
329	reflections.applyreflectionfromtheright(ref a, taubuf[i], ref t, blockstart+blocksize, m-1, blockstart+i, n-1, ref work);
330	}
331	}
332	}
333
334	//
335	// Advance
336	//
337	blockstart = blockstart+blocksize;
338	}
339	}
340
341
342	/*************************************************************************
343	QR decomposition of a rectangular complex matrix of size MxN
344
345	Input parameters:
346	A - matrix A whose indexes range within [0..M-1, 0..N-1]
347	M - number of rows in matrix A.
348	N - number of columns in matrix A.
349
350	Output parameters:
351	A - matrices Q and R in compact form
352	Tau - array of scalar factors which are used to form matrix Q. Array
353	whose indexes range within [0.. Min(M,N)-1]
354
355	Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
356	MxM, R - upper triangular (or upper trapezoid) matrix of size MxN.
357
358	-- LAPACK routine (version 3.0) --
359	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
360	Courant Institute, Argonne National Lab, and Rice University
361	September 30, 1994
362	*************************************************************************/
363	public static void cmatrixqr(ref AP.Complex[,] a,
364	int m,
365	int n,
366	ref AP.Complex[] tau)
367	{
368	AP.Complex[] work = new AP.Complex[0];
369	AP.Complex[] t = new AP.Complex[0];
370	AP.Complex[] taubuf = new AP.Complex[0];
371	int minmn = 0;
372	AP.Complex[,] tmpa = new AP.Complex[0,0];
373	AP.Complex[,] tmpt = new AP.Complex[0,0];
374	AP.Complex[,] tmpr = new AP.Complex[0,0];
375	int blockstart = 0;
376	int blocksize = 0;
377	int rowscount = 0;
378	int i = 0;
379	int j = 0;
380	int k = 0;
381	AP.Complex v = 0;
382	int i_ = 0;
383	int i1_ = 0;
384
385	if( m<=0 \| n<=0 )
386	{
387	return;
388	}
389	minmn = Math.Min(m, n);
390	work = new AP.Complex[Math.Max(m, n)+1];
391	t = new AP.Complex[Math.Max(m, n)+1];
392	tau = new AP.Complex[minmn];
393	taubuf = new AP.Complex[minmn];
394	tmpa = new AP.Complex[m, ablas.ablascomplexblocksize(ref a)];
395	tmpt = new AP.Complex[ablas.ablascomplexblocksize(ref a), ablas.ablascomplexblocksize(ref a)];
396	tmpr = new AP.Complex[2*ablas.ablascomplexblocksize(ref a), n];
397
398	//
399	// Blocked code
400	//
401	blockstart = 0;
402	while( blockstart!=minmn )
403	{
404
405	//
406	// Determine block size
407	//
408	blocksize = minmn-blockstart;
409	if( blocksize>ablas.ablascomplexblocksize(ref a) )
410	{
411	blocksize = ablas.ablascomplexblocksize(ref a);
412	}
413	rowscount = m-blockstart;
414
415	//
416	// QR decomposition of submatrix.
417	// Matrix is copied to temporary storage to solve
418	// some TLB issues arising from non-contiguous memory
419	// access pattern.
420	//
421	ablas.cmatrixcopy(rowscount, blocksize, ref a, blockstart, blockstart, ref tmpa, 0, 0);
422	cmatrixqrbasecase(ref tmpa, rowscount, blocksize, ref work, ref t, ref taubuf);
423	ablas.cmatrixcopy(rowscount, blocksize, ref tmpa, 0, 0, ref a, blockstart, blockstart);
424	i1_ = (0) - (blockstart);
425	for(i_=blockstart; i_<=blockstart+blocksize-1;i_++)
426	{
427	tau[i_] = taubuf[i_+i1_];
428	}
429
430	//
431	// Update the rest, choose between:
432	// a) Level 2 algorithm (when the rest of the matrix is small enough)
433	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
434	// representation for products of Householder transformations',
435	// by R. Schreiber and C. Van Loan.
436	//
437	if( blockstart+blocksize<=n-1 )
438	{
439	if( n-blockstart-blocksize>=2*ablas.ablascomplexblocksize(ref a) )
440	{
441
442	//
443	// Prepare block reflector
444	//
445	cmatrixblockreflector(ref tmpa, ref taubuf, true, rowscount, blocksize, ref tmpt, ref work);
446
447	//
448	// Multiply the rest of A by Q'.
449	//
450	// Q = E + YTY' = E + TmpATmpTTmpA'
451	// Q' = E + YT'Y' = E + TmpATmpT'TmpA'
452	//
453	ablas.cmatrixgemm(blocksize, n-blockstart-blocksize, rowscount, 1.0, ref tmpa, 0, 0, 2, ref a, blockstart, blockstart+blocksize, 0, 0.0, ref tmpr, 0, 0);
454	ablas.cmatrixgemm(blocksize, n-blockstart-blocksize, blocksize, 1.0, ref tmpt, 0, 0, 2, ref tmpr, 0, 0, 0, 0.0, ref tmpr, blocksize, 0);
455	ablas.cmatrixgemm(rowscount, n-blockstart-blocksize, blocksize, 1.0, ref tmpa, 0, 0, 0, ref tmpr, blocksize, 0, 0, 1.0, ref a, blockstart, blockstart+blocksize);
456	}
457	else
458	{
459
460	//
461	// Level 2 algorithm
462	//
463	for(i=0; i<=blocksize-1; i++)
464	{
465	i1_ = (i) - (1);
466	for(i_=1; i_<=rowscount-i;i_++)
467	{
468	t[i_] = tmpa[i_+i1_,i];
469	}
470	t[1] = 1;
471	creflections.complexapplyreflectionfromtheleft(ref a, AP.Math.Conj(taubuf[i]), ref t, blockstart+i, m-1, blockstart+blocksize, n-1, ref work);
472	}
473	}
474	}
475
476	//
477	// Advance
478	//
479	blockstart = blockstart+blocksize;
480	}
481	}
482
483
484	/*************************************************************************
485	LQ decomposition of a rectangular complex matrix of size MxN
486
487	Input parameters:
488	A - matrix A whose indexes range within [0..M-1, 0..N-1]
489	M - number of rows in matrix A.
490	N - number of columns in matrix A.
491
492	Output parameters:
493	A - matrices Q and L in compact form
494	Tau - array of scalar factors which are used to form matrix Q. Array
495	whose indexes range within [0.. Min(M,N)-1]
496
497	Matrix A is represented as A = LQ, where Q is an orthogonal matrix of size
498	MxM, L - lower triangular (or lower trapezoid) matrix of size MxN.
499
500	-- LAPACK routine (version 3.0) --
501	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
502	Courant Institute, Argonne National Lab, and Rice University
503	September 30, 1994
504	*************************************************************************/
505	public static void cmatrixlq(ref AP.Complex[,] a,
506	int m,
507	int n,
508	ref AP.Complex[] tau)
509	{
510	AP.Complex[] work = new AP.Complex[0];
511	AP.Complex[] t = new AP.Complex[0];
512	AP.Complex[] taubuf = new AP.Complex[0];
513	int minmn = 0;
514	AP.Complex[,] tmpa = new AP.Complex[0,0];
515	AP.Complex[,] tmpt = new AP.Complex[0,0];
516	AP.Complex[,] tmpr = new AP.Complex[0,0];
517	int blockstart = 0;
518	int blocksize = 0;
519	int columnscount = 0;
520	int i = 0;
521	int j = 0;
522	int k = 0;
523	AP.Complex v = 0;
524	int i_ = 0;
525	int i1_ = 0;
526
527	if( m<=0 \| n<=0 )
528	{
529	return;
530	}
531	minmn = Math.Min(m, n);
532	work = new AP.Complex[Math.Max(m, n)+1];
533	t = new AP.Complex[Math.Max(m, n)+1];
534	tau = new AP.Complex[minmn];
535	taubuf = new AP.Complex[minmn];
536	tmpa = new AP.Complex[ablas.ablascomplexblocksize(ref a), n];
537	tmpt = new AP.Complex[ablas.ablascomplexblocksize(ref a), ablas.ablascomplexblocksize(ref a)];
538	tmpr = new AP.Complex[m, 2*ablas.ablascomplexblocksize(ref a)];
539
540	//
541	// Blocked code
542	//
543	blockstart = 0;
544	while( blockstart!=minmn )
545	{
546
547	//
548	// Determine block size
549	//
550	blocksize = minmn-blockstart;
551	if( blocksize>ablas.ablascomplexblocksize(ref a) )
552	{
553	blocksize = ablas.ablascomplexblocksize(ref a);
554	}
555	columnscount = n-blockstart;
556
557	//
558	// LQ decomposition of submatrix.
559	// Matrix is copied to temporary storage to solve
560	// some TLB issues arising from non-contiguous memory
561	// access pattern.
562	//
563	ablas.cmatrixcopy(blocksize, columnscount, ref a, blockstart, blockstart, ref tmpa, 0, 0);
564	cmatrixlqbasecase(ref tmpa, blocksize, columnscount, ref work, ref t, ref taubuf);
565	ablas.cmatrixcopy(blocksize, columnscount, ref tmpa, 0, 0, ref a, blockstart, blockstart);
566	i1_ = (0) - (blockstart);
567	for(i_=blockstart; i_<=blockstart+blocksize-1;i_++)
568	{
569	tau[i_] = taubuf[i_+i1_];
570	}
571
572	//
573	// Update the rest, choose between:
574	// a) Level 2 algorithm (when the rest of the matrix is small enough)
575	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
576	// representation for products of Householder transformations',
577	// by R. Schreiber and C. Van Loan.
578	//
579	if( blockstart+blocksize<=m-1 )
580	{
581	if( m-blockstart-blocksize>=2*ablas.ablascomplexblocksize(ref a) )
582	{
583
584	//
585	// Prepare block reflector
586	//
587	cmatrixblockreflector(ref tmpa, ref taubuf, false, columnscount, blocksize, ref tmpt, ref work);
588
589	//
590	// Multiply the rest of A by Q.
591	//
592	// Q = E + YTY' = E + TmpA'TmpTTmpA
593	//
594	ablas.cmatrixgemm(m-blockstart-blocksize, blocksize, columnscount, 1.0, ref a, blockstart+blocksize, blockstart, 0, ref tmpa, 0, 0, 2, 0.0, ref tmpr, 0, 0);
595	ablas.cmatrixgemm(m-blockstart-blocksize, blocksize, blocksize, 1.0, ref tmpr, 0, 0, 0, ref tmpt, 0, 0, 0, 0.0, ref tmpr, 0, blocksize);
596	ablas.cmatrixgemm(m-blockstart-blocksize, columnscount, blocksize, 1.0, ref tmpr, 0, blocksize, 0, ref tmpa, 0, 0, 0, 1.0, ref a, blockstart+blocksize, blockstart);
597	}
598	else
599	{
600
601	//
602	// Level 2 algorithm
603	//
604	for(i=0; i<=blocksize-1; i++)
605	{
606	i1_ = (i) - (1);
607	for(i_=1; i_<=columnscount-i;i_++)
608	{
609	t[i_] = AP.Math.Conj(tmpa[i,i_+i1_]);
610	}
611	t[1] = 1;
612	creflections.complexapplyreflectionfromtheright(ref a, taubuf[i], ref t, blockstart+blocksize, m-1, blockstart+i, n-1, ref work);
613	}
614	}
615	}
616
617	//
618	// Advance
619	//
620	blockstart = blockstart+blocksize;
621	}
622	}
623
624
625	/*************************************************************************
626	Partial unpacking of matrix Q from the QR decomposition of a matrix A
627
628	Input parameters:
629	A - matrices Q and R in compact form.
630	Output of RMatrixQR subroutine.
631	M - number of rows in given matrix A. M>=0.
632	N - number of columns in given matrix A. N>=0.
633	Tau - scalar factors which are used to form Q.
634	Output of the RMatrixQR subroutine.
635	QColumns - required number of columns of matrix Q. M>=QColumns>=0.
636
637	Output parameters:
638	Q - first QColumns columns of matrix Q.
639	Array whose indexes range within [0..M-1, 0..QColumns-1].
640	If QColumns=0, the array remains unchanged.
641
642	-- ALGLIB routine --
643	17.02.2010
644	Bochkanov Sergey
645	*************************************************************************/
646	public static void rmatrixqrunpackq(ref double[,] a,
647	int m,
648	int n,
649	ref double[] tau,
650	int qcolumns,
651	ref double[,] q)
652	{
653	double[] work = new double[0];
654	double[] t = new double[0];
655	double[] taubuf = new double[0];
656	int minmn = 0;
657	int refcnt = 0;
658	double[,] tmpa = new double[0,0];
659	double[,] tmpt = new double[0,0];
660	double[,] tmpr = new double[0,0];
661	int blockstart = 0;
662	int blocksize = 0;
663	int rowscount = 0;
664	int i = 0;
665	int j = 0;
666	int k = 0;
667	double v = 0;
668	int i_ = 0;
669	int i1_ = 0;
670
671	System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
672	if( m<=0 \| n<=0 \| qcolumns<=0 )
673	{
674	return;
675	}
676
677	//
678	// init
679	//
680	minmn = Math.Min(m, n);
681	refcnt = Math.Min(minmn, qcolumns);
682	q = new double[m, qcolumns];
683	for(i=0; i<=m-1; i++)
684	{
685	for(j=0; j<=qcolumns-1; j++)
686	{
687	if( i==j )
688	{
689	q[i,j] = 1;
690	}
691	else
692	{
693	q[i,j] = 0;
694	}
695	}
696	}
697	work = new double[Math.Max(m, qcolumns)+1];
698	t = new double[Math.Max(m, qcolumns)+1];
699	taubuf = new double[minmn];
700	tmpa = new double[m, ablas.ablasblocksize(ref a)];
701	tmpt = new double[ablas.ablasblocksize(ref a), 2*ablas.ablasblocksize(ref a)];
702	tmpr = new double[2*ablas.ablasblocksize(ref a), qcolumns];
703
704	//
705	// Blocked code
706	//
707	blockstart = ablas.ablasblocksize(ref a)*(refcnt/ablas.ablasblocksize(ref a));
708	blocksize = refcnt-blockstart;
709	while( blockstart>=0 )
710	{
711	rowscount = m-blockstart;
712
713	//
714	// Copy current block
715	//
716	ablas.rmatrixcopy(rowscount, blocksize, ref a, blockstart, blockstart, ref tmpa, 0, 0);
717	i1_ = (blockstart) - (0);
718	for(i_=0; i_<=blocksize-1;i_++)
719	{
720	taubuf[i_] = tau[i_+i1_];
721	}
722
723	//
724	// Update, choose between:
725	// a) Level 2 algorithm (when the rest of the matrix is small enough)
726	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
727	// representation for products of Householder transformations',
728	// by R. Schreiber and C. Van Loan.
729	//
730	if( qcolumns>=2*ablas.ablasblocksize(ref a) )
731	{
732
733	//
734	// Prepare block reflector
735	//
736	rmatrixblockreflector(ref tmpa, ref taubuf, true, rowscount, blocksize, ref tmpt, ref work);
737
738	//
739	// Multiply matrix by Q.
740	//
741	// Q = E + YTY' = E + TmpATmpTTmpA'
742	//
743	ablas.rmatrixgemm(blocksize, qcolumns, rowscount, 1.0, ref tmpa, 0, 0, 1, ref q, blockstart, 0, 0, 0.0, ref tmpr, 0, 0);
744	ablas.rmatrixgemm(blocksize, qcolumns, blocksize, 1.0, ref tmpt, 0, 0, 0, ref tmpr, 0, 0, 0, 0.0, ref tmpr, blocksize, 0);
745	ablas.rmatrixgemm(rowscount, qcolumns, blocksize, 1.0, ref tmpa, 0, 0, 0, ref tmpr, blocksize, 0, 0, 1.0, ref q, blockstart, 0);
746	}
747	else
748	{
749
750	//
751	// Level 2 algorithm
752	//
753	for(i=blocksize-1; i>=0; i--)
754	{
755	i1_ = (i) - (1);
756	for(i_=1; i_<=rowscount-i;i_++)
757	{
758	t[i_] = tmpa[i_+i1_,i];
759	}
760	t[1] = 1;
761	reflections.applyreflectionfromtheleft(ref q, taubuf[i], ref t, blockstart+i, m-1, 0, qcolumns-1, ref work);
762	}
763	}
764
765	//
766	// Advance
767	//
768	blockstart = blockstart-ablas.ablasblocksize(ref a);
769	blocksize = ablas.ablasblocksize(ref a);
770	}
771	}
772
773
774	/*************************************************************************
775	Unpacking of matrix R from the QR decomposition of a matrix A
776
777	Input parameters:
778	A - matrices Q and R in compact form.
779	Output of RMatrixQR subroutine.
780	M - number of rows in given matrix A. M>=0.
781	N - number of columns in given matrix A. N>=0.
782
783	Output parameters:
784	R - matrix R, array[0..M-1, 0..N-1].
785
786	-- ALGLIB routine --
787	17.02.2010
788	Bochkanov Sergey
789	*************************************************************************/
790	public static void rmatrixqrunpackr(ref double[,] a,
791	int m,
792	int n,
793	ref double[,] r)
794	{
795	int i = 0;
796	int k = 0;
797	int i_ = 0;
798
799	if( m<=0 \| n<=0 )
800	{
801	return;
802	}
803	k = Math.Min(m, n);
804	r = new double[m, n];
805	for(i=0; i<=n-1; i++)
806	{
807	r[0,i] = 0;
808	}
809	for(i=1; i<=m-1; i++)
810	{
811	for(i_=0; i_<=n-1;i_++)
812	{
813	r[i,i_] = r[0,i_];
814	}
815	}
816	for(i=0; i<=k-1; i++)
817	{
818	for(i_=i; i_<=n-1;i_++)
819	{
820	r[i,i_] = a[i,i_];
821	}
822	}
823	}
824
825
826	/*************************************************************************
827	Partial unpacking of matrix Q from the LQ decomposition of a matrix A
828
829	Input parameters:
830	A - matrices L and Q in compact form.
831	Output of RMatrixLQ subroutine.
832	M - number of rows in given matrix A. M>=0.
833	N - number of columns in given matrix A. N>=0.
834	Tau - scalar factors which are used to form Q.
835	Output of the RMatrixLQ subroutine.
836	QRows - required number of rows in matrix Q. N>=QRows>=0.
837
838	Output parameters:
839	Q - first QRows rows of matrix Q. Array whose indexes range
840	within [0..QRows-1, 0..N-1]. If QRows=0, the array remains
841	unchanged.
842
843	-- ALGLIB routine --
844	17.02.2010
845	Bochkanov Sergey
846	*************************************************************************/
847	public static void rmatrixlqunpackq(ref double[,] a,
848	int m,
849	int n,
850	ref double[] tau,
851	int qrows,
852	ref double[,] q)
853	{
854	double[] work = new double[0];
855	double[] t = new double[0];
856	double[] taubuf = new double[0];
857	int minmn = 0;
858	int refcnt = 0;
859	double[,] tmpa = new double[0,0];
860	double[,] tmpt = new double[0,0];
861	double[,] tmpr = new double[0,0];
862	int blockstart = 0;
863	int blocksize = 0;
864	int columnscount = 0;
865	int i = 0;
866	int j = 0;
867	int k = 0;
868	double v = 0;
869	int i_ = 0;
870	int i1_ = 0;
871
872	System.Diagnostics.Debug.Assert(qrows<=n, "RMatrixLQUnpackQ: QRows>N!");
873	if( m<=0 \| n<=0 \| qrows<=0 )
874	{
875	return;
876	}
877
878	//
879	// init
880	//
881	minmn = Math.Min(m, n);
882	refcnt = Math.Min(minmn, qrows);
883	work = new double[Math.Max(m, n)+1];
884	t = new double[Math.Max(m, n)+1];
885	taubuf = new double[minmn];
886	tmpa = new double[ablas.ablasblocksize(ref a), n];
887	tmpt = new double[ablas.ablasblocksize(ref a), 2*ablas.ablasblocksize(ref a)];
888	tmpr = new double[qrows, 2*ablas.ablasblocksize(ref a)];
889	q = new double[qrows, n];
890	for(i=0; i<=qrows-1; i++)
891	{
892	for(j=0; j<=n-1; j++)
893	{
894	if( i==j )
895	{
896	q[i,j] = 1;
897	}
898	else
899	{
900	q[i,j] = 0;
901	}
902	}
903	}
904
905	//
906	// Blocked code
907	//
908	blockstart = ablas.ablasblocksize(ref a)*(refcnt/ablas.ablasblocksize(ref a));
909	blocksize = refcnt-blockstart;
910	while( blockstart>=0 )
911	{
912	columnscount = n-blockstart;
913
914	//
915	// Copy submatrix
916	//
917	ablas.rmatrixcopy(blocksize, columnscount, ref a, blockstart, blockstart, ref tmpa, 0, 0);
918	i1_ = (blockstart) - (0);
919	for(i_=0; i_<=blocksize-1;i_++)
920	{
921	taubuf[i_] = tau[i_+i1_];
922	}
923
924	//
925	// Update matrix, choose between:
926	// a) Level 2 algorithm (when the rest of the matrix is small enough)
927	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
928	// representation for products of Householder transformations',
929	// by R. Schreiber and C. Van Loan.
930	//
931	if( qrows>=2*ablas.ablasblocksize(ref a) )
932	{
933
934	//
935	// Prepare block reflector
936	//
937	rmatrixblockreflector(ref tmpa, ref taubuf, false, columnscount, blocksize, ref tmpt, ref work);
938
939	//
940	// Multiply the rest of A by Q'.
941	//
942	// Q' = E + YT'Y' = E + TmpA'TmpT'TmpA
943	//
944	ablas.rmatrixgemm(qrows, blocksize, columnscount, 1.0, ref q, 0, blockstart, 0, ref tmpa, 0, 0, 1, 0.0, ref tmpr, 0, 0);
945	ablas.rmatrixgemm(qrows, blocksize, blocksize, 1.0, ref tmpr, 0, 0, 0, ref tmpt, 0, 0, 1, 0.0, ref tmpr, 0, blocksize);
946	ablas.rmatrixgemm(qrows, columnscount, blocksize, 1.0, ref tmpr, 0, blocksize, 0, ref tmpa, 0, 0, 0, 1.0, ref q, 0, blockstart);
947	}
948	else
949	{
950
951	//
952	// Level 2 algorithm
953	//
954	for(i=blocksize-1; i>=0; i--)
955	{
956	i1_ = (i) - (1);
957	for(i_=1; i_<=columnscount-i;i_++)
958	{
959	t[i_] = tmpa[i,i_+i1_];
960	}
961	t[1] = 1;
962	reflections.applyreflectionfromtheright(ref q, taubuf[i], ref t, 0, qrows-1, blockstart+i, n-1, ref work);
963	}
964	}
965
966	//
967	// Advance
968	//
969	blockstart = blockstart-ablas.ablasblocksize(ref a);
970	blocksize = ablas.ablasblocksize(ref a);
971	}
972	}
973
974
975	/*************************************************************************
976	Unpacking of matrix L from the LQ decomposition of a matrix A
977
978	Input parameters:
979	A - matrices Q and L in compact form.
980	Output of RMatrixLQ subroutine.
981	M - number of rows in given matrix A. M>=0.
982	N - number of columns in given matrix A. N>=0.
983
984	Output parameters:
985	L - matrix L, array[0..M-1, 0..N-1].
986
987	-- ALGLIB routine --
988	17.02.2010
989	Bochkanov Sergey
990	*************************************************************************/
991	public static void rmatrixlqunpackl(ref double[,] a,
992	int m,
993	int n,
994	ref double[,] l)
995	{
996	int i = 0;
997	int k = 0;
998	int i_ = 0;
999
1000	if( m<=0 \| n<=0 )
1001	{
1002	return;
1003	}
1004	l = new double[m, n];
1005	for(i=0; i<=n-1; i++)
1006	{
1007	l[0,i] = 0;
1008	}
1009	for(i=1; i<=m-1; i++)
1010	{
1011	for(i_=0; i_<=n-1;i_++)
1012	{
1013	l[i,i_] = l[0,i_];
1014	}
1015	}
1016	for(i=0; i<=m-1; i++)
1017	{
1018	k = Math.Min(i, n-1);
1019	for(i_=0; i_<=k;i_++)
1020	{
1021	l[i,i_] = a[i,i_];
1022	}
1023	}
1024	}
1025
1026
1027	/*************************************************************************
1028	Partial unpacking of matrix Q from QR decomposition of a complex matrix A.
1029
1030	Input parameters:
1031	A - matrices Q and R in compact form.
1032	Output of CMatrixQR subroutine .
1033	M - number of rows in matrix A. M>=0.
1034	N - number of columns in matrix A. N>=0.
1035	Tau - scalar factors which are used to form Q.
1036	Output of CMatrixQR subroutine .
1037	QColumns - required number of columns in matrix Q. M>=QColumns>=0.
1038
1039	Output parameters:
1040	Q - first QColumns columns of matrix Q.
1041	Array whose index ranges within [0..M-1, 0..QColumns-1].
1042	If QColumns=0, array isn't changed.
1043
1044	-- ALGLIB routine --
1045	17.02.2010
1046	Bochkanov Sergey
1047	*************************************************************************/
1048	public static void cmatrixqrunpackq(ref AP.Complex[,] a,
1049	int m,
1050	int n,
1051	ref AP.Complex[] tau,
1052	int qcolumns,
1053	ref AP.Complex[,] q)
1054	{
1055	AP.Complex[] work = new AP.Complex[0];
1056	AP.Complex[] t = new AP.Complex[0];
1057	AP.Complex[] taubuf = new AP.Complex[0];
1058	int minmn = 0;
1059	int refcnt = 0;
1060	AP.Complex[,] tmpa = new AP.Complex[0,0];
1061	AP.Complex[,] tmpt = new AP.Complex[0,0];
1062	AP.Complex[,] tmpr = new AP.Complex[0,0];
1063	int blockstart = 0;
1064	int blocksize = 0;
1065	int rowscount = 0;
1066	int i = 0;
1067	int j = 0;
1068	int k = 0;
1069	AP.Complex v = 0;
1070	int i_ = 0;
1071	int i1_ = 0;
1072
1073	System.Diagnostics.Debug.Assert(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
1074	if( m<=0 \| n<=0 )
1075	{
1076	return;
1077	}
1078
1079	//
1080	// init
1081	//
1082	minmn = Math.Min(m, n);
1083	refcnt = Math.Min(minmn, qcolumns);
1084	work = new AP.Complex[Math.Max(m, n)+1];
1085	t = new AP.Complex[Math.Max(m, n)+1];
1086	taubuf = new AP.Complex[minmn];
1087	tmpa = new AP.Complex[m, ablas.ablascomplexblocksize(ref a)];
1088	tmpt = new AP.Complex[ablas.ablascomplexblocksize(ref a), ablas.ablascomplexblocksize(ref a)];
1089	tmpr = new AP.Complex[2*ablas.ablascomplexblocksize(ref a), qcolumns];
1090	q = new AP.Complex[m, qcolumns];
1091	for(i=0; i<=m-1; i++)
1092	{
1093	for(j=0; j<=qcolumns-1; j++)
1094	{
1095	if( i==j )
1096	{
1097	q[i,j] = 1;
1098	}
1099	else
1100	{
1101	q[i,j] = 0;
1102	}
1103	}
1104	}
1105
1106	//
1107	// Blocked code
1108	//
1109	blockstart = ablas.ablascomplexblocksize(ref a)*(refcnt/ablas.ablascomplexblocksize(ref a));
1110	blocksize = refcnt-blockstart;
1111	while( blockstart>=0 )
1112	{
1113	rowscount = m-blockstart;
1114
1115	//
1116	// QR decomposition of submatrix.
1117	// Matrix is copied to temporary storage to solve
1118	// some TLB issues arising from non-contiguous memory
1119	// access pattern.
1120	//
1121	ablas.cmatrixcopy(rowscount, blocksize, ref a, blockstart, blockstart, ref tmpa, 0, 0);
1122	i1_ = (blockstart) - (0);
1123	for(i_=0; i_<=blocksize-1;i_++)
1124	{
1125	taubuf[i_] = tau[i_+i1_];
1126	}
1127
1128	//
1129	// Update matrix, choose between:
1130	// a) Level 2 algorithm (when the rest of the matrix is small enough)
1131	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
1132	// representation for products of Householder transformations',
1133	// by R. Schreiber and C. Van Loan.
1134	//
1135	if( qcolumns>=2*ablas.ablascomplexblocksize(ref a) )
1136	{
1137
1138	//
1139	// Prepare block reflector
1140	//
1141	cmatrixblockreflector(ref tmpa, ref taubuf, true, rowscount, blocksize, ref tmpt, ref work);
1142
1143	//
1144	// Multiply the rest of A by Q.
1145	//
1146	// Q = E + YTY' = E + TmpATmpTTmpA'
1147	//
1148	ablas.cmatrixgemm(blocksize, qcolumns, rowscount, 1.0, ref tmpa, 0, 0, 2, ref q, blockstart, 0, 0, 0.0, ref tmpr, 0, 0);
1149	ablas.cmatrixgemm(blocksize, qcolumns, blocksize, 1.0, ref tmpt, 0, 0, 0, ref tmpr, 0, 0, 0, 0.0, ref tmpr, blocksize, 0);
1150	ablas.cmatrixgemm(rowscount, qcolumns, blocksize, 1.0, ref tmpa, 0, 0, 0, ref tmpr, blocksize, 0, 0, 1.0, ref q, blockstart, 0);
1151	}
1152	else
1153	{
1154
1155	//
1156	// Level 2 algorithm
1157	//
1158	for(i=blocksize-1; i>=0; i--)
1159	{
1160	i1_ = (i) - (1);
1161	for(i_=1; i_<=rowscount-i;i_++)
1162	{
1163	t[i_] = tmpa[i_+i1_,i];
1164	}
1165	t[1] = 1;
1166	creflections.complexapplyreflectionfromtheleft(ref q, taubuf[i], ref t, blockstart+i, m-1, 0, qcolumns-1, ref work);
1167	}
1168	}
1169
1170	//
1171	// Advance
1172	//
1173	blockstart = blockstart-ablas.ablascomplexblocksize(ref a);
1174	blocksize = ablas.ablascomplexblocksize(ref a);
1175	}
1176	}
1177
1178
1179	/*************************************************************************
1180	Unpacking of matrix R from the QR decomposition of a matrix A
1181
1182	Input parameters:
1183	A - matrices Q and R in compact form.
1184	Output of CMatrixQR subroutine.
1185	M - number of rows in given matrix A. M>=0.
1186	N - number of columns in given matrix A. N>=0.
1187
1188	Output parameters:
1189	R - matrix R, array[0..M-1, 0..N-1].
1190
1191	-- ALGLIB routine --
1192	17.02.2010
1193	Bochkanov Sergey
1194	*************************************************************************/
1195	public static void cmatrixqrunpackr(ref AP.Complex[,] a,
1196	int m,
1197	int n,
1198	ref AP.Complex[,] r)
1199	{
1200	int i = 0;
1201	int k = 0;
1202	int i_ = 0;
1203
1204	if( m<=0 \| n<=0 )
1205	{
1206	return;
1207	}
1208	k = Math.Min(m, n);
1209	r = new AP.Complex[m, n];
1210	for(i=0; i<=n-1; i++)
1211	{
1212	r[0,i] = 0;
1213	}
1214	for(i=1; i<=m-1; i++)
1215	{
1216	for(i_=0; i_<=n-1;i_++)
1217	{
1218	r[i,i_] = r[0,i_];
1219	}
1220	}
1221	for(i=0; i<=k-1; i++)
1222	{
1223	for(i_=i; i_<=n-1;i_++)
1224	{
1225	r[i,i_] = a[i,i_];
1226	}
1227	}
1228	}
1229
1230
1231	/*************************************************************************
1232	Partial unpacking of matrix Q from LQ decomposition of a complex matrix A.
1233
1234	Input parameters:
1235	A - matrices Q and R in compact form.
1236	Output of CMatrixLQ subroutine .
1237	M - number of rows in matrix A. M>=0.
1238	N - number of columns in matrix A. N>=0.
1239	Tau - scalar factors which are used to form Q.
1240	Output of CMatrixLQ subroutine .
1241	QRows - required number of rows in matrix Q. N>=QColumns>=0.
1242
1243	Output parameters:
1244	Q - first QRows rows of matrix Q.
1245	Array whose index ranges within [0..QRows-1, 0..N-1].
1246	If QRows=0, array isn't changed.
1247
1248	-- ALGLIB routine --
1249	17.02.2010
1250	Bochkanov Sergey
1251	*************************************************************************/
1252	public static void cmatrixlqunpackq(ref AP.Complex[,] a,
1253	int m,
1254	int n,
1255	ref AP.Complex[] tau,
1256	int qrows,
1257	ref AP.Complex[,] q)
1258	{
1259	AP.Complex[] work = new AP.Complex[0];
1260	AP.Complex[] t = new AP.Complex[0];
1261	AP.Complex[] taubuf = new AP.Complex[0];
1262	int minmn = 0;
1263	int refcnt = 0;
1264	AP.Complex[,] tmpa = new AP.Complex[0,0];
1265	AP.Complex[,] tmpt = new AP.Complex[0,0];
1266	AP.Complex[,] tmpr = new AP.Complex[0,0];
1267	int blockstart = 0;
1268	int blocksize = 0;
1269	int columnscount = 0;
1270	int i = 0;
1271	int j = 0;
1272	int k = 0;
1273	AP.Complex v = 0;
1274	int i_ = 0;
1275	int i1_ = 0;
1276
1277	if( m<=0 \| n<=0 )
1278	{
1279	return;
1280	}
1281
1282	//
1283	// Init
1284	//
1285	minmn = Math.Min(m, n);
1286	refcnt = Math.Min(minmn, qrows);
1287	work = new AP.Complex[Math.Max(m, n)+1];
1288	t = new AP.Complex[Math.Max(m, n)+1];
1289	taubuf = new AP.Complex[minmn];
1290	tmpa = new AP.Complex[ablas.ablascomplexblocksize(ref a), n];
1291	tmpt = new AP.Complex[ablas.ablascomplexblocksize(ref a), ablas.ablascomplexblocksize(ref a)];
1292	tmpr = new AP.Complex[qrows, 2*ablas.ablascomplexblocksize(ref a)];
1293	q = new AP.Complex[qrows, n];
1294	for(i=0; i<=qrows-1; i++)
1295	{
1296	for(j=0; j<=n-1; j++)
1297	{
1298	if( i==j )
1299	{
1300	q[i,j] = 1;
1301	}
1302	else
1303	{
1304	q[i,j] = 0;
1305	}
1306	}
1307	}
1308
1309	//
1310	// Blocked code
1311	//
1312	blockstart = ablas.ablascomplexblocksize(ref a)*(refcnt/ablas.ablascomplexblocksize(ref a));
1313	blocksize = refcnt-blockstart;
1314	while( blockstart>=0 )
1315	{
1316	columnscount = n-blockstart;
1317
1318	//
1319	// LQ decomposition of submatrix.
1320	// Matrix is copied to temporary storage to solve
1321	// some TLB issues arising from non-contiguous memory
1322	// access pattern.
1323	//
1324	ablas.cmatrixcopy(blocksize, columnscount, ref a, blockstart, blockstart, ref tmpa, 0, 0);
1325	i1_ = (blockstart) - (0);
1326	for(i_=0; i_<=blocksize-1;i_++)
1327	{
1328	taubuf[i_] = tau[i_+i1_];
1329	}
1330
1331	//
1332	// Update matrix, choose between:
1333	// a) Level 2 algorithm (when the rest of the matrix is small enough)
1334	// b) blocked algorithm, see algorithm 5 from 'A storage efficient WY
1335	// representation for products of Householder transformations',
1336	// by R. Schreiber and C. Van Loan.
1337	//
1338	if( qrows>=2*ablas.ablascomplexblocksize(ref a) )
1339	{
1340
1341	//
1342	// Prepare block reflector
1343	//
1344	cmatrixblockreflector(ref tmpa, ref taubuf, false, columnscount, blocksize, ref tmpt, ref work);
1345
1346	//
1347	// Multiply the rest of A by Q'.
1348	//
1349	// Q' = E + YT'Y' = E + TmpA'TmpT'TmpA
1350	//
1351	ablas.cmatrixgemm(qrows, blocksize, columnscount, 1.0, ref q, 0, blockstart, 0, ref tmpa, 0, 0, 2, 0.0, ref tmpr, 0, 0);
1352	ablas.cmatrixgemm(qrows, blocksize, blocksize, 1.0, ref tmpr, 0, 0, 0, ref tmpt, 0, 0, 2, 0.0, ref tmpr, 0, blocksize);
1353	ablas.cmatrixgemm(qrows, columnscount, blocksize, 1.0, ref tmpr, 0, blocksize, 0, ref tmpa, 0, 0, 0, 1.0, ref q, 0, blockstart);
1354	}
1355	else
1356	{
1357
1358	//
1359	// Level 2 algorithm
1360	//
1361	for(i=blocksize-1; i>=0; i--)
1362	{
1363	i1_ = (i) - (1);
1364	for(i_=1; i_<=columnscount-i;i_++)
1365	{
1366	t[i_] = AP.Math.Conj(tmpa[i,i_+i1_]);
1367	}
1368	t[1] = 1;
1369	creflections.complexapplyreflectionfromtheright(ref q, AP.Math.Conj(taubuf[i]), ref t, 0, qrows-1, blockstart+i, n-1, ref work);
1370	}
1371	}
1372
1373	//
1374	// Advance
1375	//
1376	blockstart = blockstart-ablas.ablascomplexblocksize(ref a);
1377	blocksize = ablas.ablascomplexblocksize(ref a);
1378	}
1379	}
1380
1381
1382	/*************************************************************************
1383	Unpacking of matrix L from the LQ decomposition of a matrix A
1384
1385	Input parameters:
1386	A - matrices Q and L in compact form.
1387	Output of CMatrixLQ subroutine.
1388	M - number of rows in given matrix A. M>=0.
1389	N - number of columns in given matrix A. N>=0.
1390
1391	Output parameters:
1392	L - matrix L, array[0..M-1, 0..N-1].
1393
1394	-- ALGLIB routine --
1395	17.02.2010
1396	Bochkanov Sergey
1397	*************************************************************************/
1398	public static void cmatrixlqunpackl(ref AP.Complex[,] a,
1399	int m,
1400	int n,
1401	ref AP.Complex[,] l)
1402	{
1403	int i = 0;
1404	int k = 0;
1405	int i_ = 0;
1406
1407	if( m<=0 \| n<=0 )
1408	{
1409	return;
1410	}
1411	l = new AP.Complex[m, n];
1412	for(i=0; i<=n-1; i++)
1413	{
1414	l[0,i] = 0;
1415	}
1416	for(i=1; i<=m-1; i++)
1417	{
1418	for(i_=0; i_<=n-1;i_++)
1419	{
1420	l[i,i_] = l[0,i_];
1421	}
1422	}
1423	for(i=0; i<=m-1; i++)
1424	{
1425	k = Math.Min(i, n-1);
1426	for(i_=0; i_<=k;i_++)
1427	{
1428	l[i,i_] = a[i,i_];
1429	}
1430	}
1431	}
1432
1433
1434	/*************************************************************************
1435	Reduction of a rectangular matrix to bidiagonal form
1436
1437	The algorithm reduces the rectangular matrix A to bidiagonal form by
1438	orthogonal transformations P and Q: A = QBP.
1439
1440	Input parameters:
1441	A - source matrix. array[0..M-1, 0..N-1]
1442	M - number of rows in matrix A.
1443	N - number of columns in matrix A.
1444
1445	Output parameters:
1446	A - matrices Q, B, P in compact form (see below).
1447	TauQ - scalar factors which are used to form matrix Q.
1448	TauP - scalar factors which are used to form matrix P.
1449
1450	The main diagonal and one of the secondary diagonals of matrix A are
1451	replaced with bidiagonal matrix B. Other elements contain elementary
1452	reflections which form MxM matrix Q and NxN matrix P, respectively.
1453
1454	If M>=N, B is the upper bidiagonal MxN matrix and is stored in the
1455	corresponding elements of matrix A. Matrix Q is represented as a
1456	product of elementary reflections Q = H(0)H(1)...*H(n-1), where
1457	H(i) = 1-tauvv'. Here tau is a scalar which is stored in TauQ[i], and
1458	vector v has the following structure: v(0:i-1)=0, v(i)=1, v(i+1:m-1) is
1459	stored in elements A(i+1:m-1,i). Matrix P is as follows: P =
1460	G(0)G(1)...G(n-2), where G(i) = 1 - tauu*u'. Tau is stored in TauP[i],
1461	u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).
1462
1463	If M<N, B is the lower bidiagonal MxN matrix and is stored in the
1464	corresponding elements of matrix A. Q = H(0)H(1)...*H(m-2), where
1465	H(i) = 1 - tauvv', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
1466	is stored in elements A(i+2:m-1,i). P = G(0)G(1)...*G(m-1),
1467	G(i) = 1-tauuu', tau is stored in TauP, u(0:i-1)=0, u(i)=1, u(i+1:n-1)
1468	is stored in A(i,i+1:n-1).
1469
1470	EXAMPLE:
1471
1472	m=6, n=5 (m > n): m=5, n=6 (m < n):
1473
1474	( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
1475	( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
1476	( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
1477	( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
1478	( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
1479	( v1 v2 v3 v4 v5 )
1480
1481	Here vi and ui are vectors which form H(i) and G(i), and d and e -
1482	are the diagonal and off-diagonal elements of matrix B.
1483
1484	-- LAPACK routine (version 3.0) --
1485	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1486	Courant Institute, Argonne National Lab, and Rice University
1487	September 30, 1994.
1488	Sergey Bochkanov, ALGLIB project, translation from FORTRAN to
1489	pseudocode, 2007-2010.
1490	*************************************************************************/
1491	public static void rmatrixbd(ref double[,] a,
1492	int m,
1493	int n,
1494	ref double[] tauq,
1495	ref double[] taup)
1496	{
1497	double[] work = new double[0];
1498	double[] t = new double[0];
1499	int minmn = 0;
1500	int maxmn = 0;
1501	int i = 0;
1502	double ltau = 0;
1503	int i_ = 0;
1504	int i1_ = 0;
1505
1506
1507	//
1508	// Prepare
1509	//
1510	if( n<=0 \| m<=0 )
1511	{
1512	return;
1513	}
1514	minmn = Math.Min(m, n);
1515	maxmn = Math.Max(m, n);
1516	work = new double[maxmn+1];
1517	t = new double[maxmn+1];
1518	if( m>=n )
1519	{
1520	tauq = new double[n];
1521	taup = new double[n];
1522	}
1523	else
1524	{
1525	tauq = new double[m];
1526	taup = new double[m];
1527	}
1528	if( m>=n )
1529	{
1530
1531	//
1532	// Reduce to upper bidiagonal form
1533	//
1534	for(i=0; i<=n-1; i++)
1535	{
1536
1537	//
1538	// Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
1539	//
1540	i1_ = (i) - (1);
1541	for(i_=1; i_<=m-i;i_++)
1542	{
1543	t[i_] = a[i_+i1_,i];
1544	}
1545	reflections.generatereflection(ref t, m-i, ref ltau);
1546	tauq[i] = ltau;
1547	i1_ = (1) - (i);
1548	for(i_=i; i_<=m-1;i_++)
1549	{
1550	a[i_,i] = t[i_+i1_];
1551	}
1552	t[1] = 1;
1553
1554	//
1555	// Apply H(i) to A(i:m-1,i+1:n-1) from the left
1556	//
1557	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i, m-1, i+1, n-1, ref work);
1558	if( i<n-1 )
1559	{
1560
1561	//
1562	// Generate elementary reflector G(i) to annihilate
1563	// A(i,i+2:n-1)
1564	//
1565	i1_ = (i+1) - (1);
1566	for(i_=1; i_<=n-i-1;i_++)
1567	{
1568	t[i_] = a[i,i_+i1_];
1569	}
1570	reflections.generatereflection(ref t, n-1-i, ref ltau);
1571	taup[i] = ltau;
1572	i1_ = (1) - (i+1);
1573	for(i_=i+1; i_<=n-1;i_++)
1574	{
1575	a[i,i_] = t[i_+i1_];
1576	}
1577	t[1] = 1;
1578
1579	//
1580	// Apply G(i) to A(i+1:m-1,i+1:n-1) from the right
1581	//
1582	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
1583	}
1584	else
1585	{
1586	taup[i] = 0;
1587	}
1588	}
1589	}
1590	else
1591	{
1592
1593	//
1594	// Reduce to lower bidiagonal form
1595	//
1596	for(i=0; i<=m-1; i++)
1597	{
1598
1599	//
1600	// Generate elementary reflector G(i) to annihilate A(i,i+1:n-1)
1601	//
1602	i1_ = (i) - (1);
1603	for(i_=1; i_<=n-i;i_++)
1604	{
1605	t[i_] = a[i,i_+i1_];
1606	}
1607	reflections.generatereflection(ref t, n-i, ref ltau);
1608	taup[i] = ltau;
1609	i1_ = (1) - (i);
1610	for(i_=i; i_<=n-1;i_++)
1611	{
1612	a[i,i_] = t[i_+i1_];
1613	}
1614	t[1] = 1;
1615
1616	//
1617	// Apply G(i) to A(i+1:m-1,i:n-1) from the right
1618	//
1619	reflections.applyreflectionfromtheright(ref a, ltau, ref t, i+1, m-1, i, n-1, ref work);
1620	if( i<m-1 )
1621	{
1622
1623	//
1624	// Generate elementary reflector H(i) to annihilate
1625	// A(i+2:m-1,i)
1626	//
1627	i1_ = (i+1) - (1);
1628	for(i_=1; i_<=m-1-i;i_++)
1629	{
1630	t[i_] = a[i_+i1_,i];
1631	}
1632	reflections.generatereflection(ref t, m-1-i, ref ltau);
1633	tauq[i] = ltau;
1634	i1_ = (1) - (i+1);
1635	for(i_=i+1; i_<=m-1;i_++)
1636	{
1637	a[i_,i] = t[i_+i1_];
1638	}
1639	t[1] = 1;
1640
1641	//
1642	// Apply H(i) to A(i+1:m-1,i+1:n-1) from the left
1643	//
1644	reflections.applyreflectionfromtheleft(ref a, ltau, ref t, i+1, m-1, i+1, n-1, ref work);
1645	}
1646	else
1647	{
1648	tauq[i] = 0;
1649	}
1650	}
1651	}
1652	}
1653
1654
1655	/*************************************************************************
1656	Unpacking matrix Q which reduces a matrix to bidiagonal form.
1657
1658	Input parameters:
1659	QP - matrices Q and P in compact form.
1660	Output of ToBidiagonal subroutine.
1661	M - number of rows in matrix A.
1662	N - number of columns in matrix A.
1663	TAUQ - scalar factors which are used to form Q.
1664	Output of ToBidiagonal subroutine.
1665	QColumns - required number of columns in matrix Q.
1666	M>=QColumns>=0.
1667
1668	Output parameters:
1669	Q - first QColumns columns of matrix Q.
1670	Array[0..M-1, 0..QColumns-1]
1671	If QColumns=0, the array is not modified.
1672
1673	-- ALGLIB --
1674	2005-2010
1675	Bochkanov Sergey
1676	*************************************************************************/
1677	public static void rmatrixbdunpackq(ref double[,] qp,
1678	int m,
1679	int n,
1680	ref double[] tauq,
1681	int qcolumns,
1682	ref double[,] q)
1683	{
1684	int i = 0;
1685	int j = 0;
1686
1687	System.Diagnostics.Debug.Assert(qcolumns<=m, "RMatrixBDUnpackQ: QColumns>M!");
1688	System.Diagnostics.Debug.Assert(qcolumns>=0, "RMatrixBDUnpackQ: QColumns<0!");
1689	if( m==0 \| n==0 \| qcolumns==0 )
1690	{
1691	return;
1692	}
1693
1694	//
1695	// prepare Q
1696	//
1697	q = new double[m, qcolumns];
1698	for(i=0; i<=m-1; i++)
1699	{
1700	for(j=0; j<=qcolumns-1; j++)
1701	{
1702	if( i==j )
1703	{
1704	q[i,j] = 1;
1705	}
1706	else
1707	{
1708	q[i,j] = 0;
1709	}
1710	}
1711	}
1712
1713	//
1714	// Calculate
1715	//
1716	rmatrixbdmultiplybyq(ref qp, m, n, ref tauq, ref q, m, qcolumns, false, false);
1717	}
1718
1719
1720	/*************************************************************************
1721	Multiplication by matrix Q which reduces matrix A to bidiagonal form.
1722
1723	The algorithm allows pre- or post-multiply by Q or Q'.
1724
1725	Input parameters:
1726	QP - matrices Q and P in compact form.
1727	Output of ToBidiagonal subroutine.
1728	M - number of rows in matrix A.
1729	N - number of columns in matrix A.
1730	TAUQ - scalar factors which are used to form Q.
1731	Output of ToBidiagonal subroutine.
1732	Z - multiplied matrix.
1733	array[0..ZRows-1,0..ZColumns-1]
1734	ZRows - number of rows in matrix Z. If FromTheRight=False,
1735	ZRows=M, otherwise ZRows can be arbitrary.
1736	ZColumns - number of columns in matrix Z. If FromTheRight=True,
1737	ZColumns=M, otherwise ZColumns can be arbitrary.
1738	FromTheRight - pre- or post-multiply.
1739	DoTranspose - multiply by Q or Q'.
1740
1741	Output parameters:
1742	Z - product of Z and Q.
1743	Array[0..ZRows-1,0..ZColumns-1]
1744	If ZRows=0 or ZColumns=0, the array is not modified.
1745
1746	-- ALGLIB --
1747	2005-2010
1748	Bochkanov Sergey
1749	*************************************************************************/
1750	public static void rmatrixbdmultiplybyq(ref double[,] qp,
1751	int m,
1752	int n,
1753	ref double[] tauq,
1754	ref double[,] z,
1755	int zrows,
1756	int zcolumns,
1757	bool fromtheright,
1758	bool dotranspose)
1759	{
1760	int i = 0;
1761	int i1 = 0;
1762	int i2 = 0;
1763	int istep = 0;
1764	double[] v = new double[0];
1765	double[] work = new double[0];
1766	int mx = 0;
1767	int i_ = 0;
1768	int i1_ = 0;
1769
1770	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
1771	{
1772	return;
1773	}
1774	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==m \| !fromtheright & zrows==m, "RMatrixBDMultiplyByQ: incorrect Z size!");
1775
1776	//
1777	// init
1778	//
1779	mx = Math.Max(m, n);
1780	mx = Math.Max(mx, zrows);
1781	mx = Math.Max(mx, zcolumns);
1782	v = new double[mx+1];
1783	work = new double[mx+1];
1784	if( m>=n )
1785	{
1786
1787	//
1788	// setup
1789	//
1790	if( fromtheright )
1791	{
1792	i1 = 0;
1793	i2 = n-1;
1794	istep = +1;
1795	}
1796	else
1797	{
1798	i1 = n-1;
1799	i2 = 0;
1800	istep = -1;
1801	}
1802	if( dotranspose )
1803	{
1804	i = i1;
1805	i1 = i2;
1806	i2 = i;
1807	istep = -istep;
1808	}
1809
1810	//
1811	// Process
1812	//
1813	i = i1;
1814	do
1815	{
1816	i1_ = (i) - (1);
1817	for(i_=1; i_<=m-i;i_++)
1818	{
1819	v[i_] = qp[i_+i1_,i];
1820	}
1821	v[1] = 1;
1822	if( fromtheright )
1823	{
1824	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i, m-1, ref work);
1825	}
1826	else
1827	{
1828	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i, m-1, 0, zcolumns-1, ref work);
1829	}
1830	i = i+istep;
1831	}
1832	while( i!=i2+istep );
1833	}
1834	else
1835	{
1836
1837	//
1838	// setup
1839	//
1840	if( fromtheright )
1841	{
1842	i1 = 0;
1843	i2 = m-2;
1844	istep = +1;
1845	}
1846	else
1847	{
1848	i1 = m-2;
1849	i2 = 0;
1850	istep = -1;
1851	}
1852	if( dotranspose )
1853	{
1854	i = i1;
1855	i1 = i2;
1856	i2 = i;
1857	istep = -istep;
1858	}
1859
1860	//
1861	// Process
1862	//
1863	if( m-1>0 )
1864	{
1865	i = i1;
1866	do
1867	{
1868	i1_ = (i+1) - (1);
1869	for(i_=1; i_<=m-i-1;i_++)
1870	{
1871	v[i_] = qp[i_+i1_,i];
1872	}
1873	v[1] = 1;
1874	if( fromtheright )
1875	{
1876	reflections.applyreflectionfromtheright(ref z, tauq[i], ref v, 0, zrows-1, i+1, m-1, ref work);
1877	}
1878	else
1879	{
1880	reflections.applyreflectionfromtheleft(ref z, tauq[i], ref v, i+1, m-1, 0, zcolumns-1, ref work);
1881	}
1882	i = i+istep;
1883	}
1884	while( i!=i2+istep );
1885	}
1886	}
1887	}
1888
1889
1890	/*************************************************************************
1891	Unpacking matrix P which reduces matrix A to bidiagonal form.
1892	The subroutine returns transposed matrix P.
1893
1894	Input parameters:
1895	QP - matrices Q and P in compact form.
1896	Output of ToBidiagonal subroutine.
1897	M - number of rows in matrix A.
1898	N - number of columns in matrix A.
1899	TAUP - scalar factors which are used to form P.
1900	Output of ToBidiagonal subroutine.
1901	PTRows - required number of rows of matrix P^T. N >= PTRows >= 0.
1902
1903	Output parameters:
1904	PT - first PTRows columns of matrix P^T
1905	Array[0..PTRows-1, 0..N-1]
1906	If PTRows=0, the array is not modified.
1907
1908	-- ALGLIB --
1909	2005-2010
1910	Bochkanov Sergey
1911	*************************************************************************/
1912	public static void rmatrixbdunpackpt(ref double[,] qp,
1913	int m,
1914	int n,
1915	ref double[] taup,
1916	int ptrows,
1917	ref double[,] pt)
1918	{
1919	int i = 0;
1920	int j = 0;
1921
1922	System.Diagnostics.Debug.Assert(ptrows<=n, "RMatrixBDUnpackPT: PTRows>N!");
1923	System.Diagnostics.Debug.Assert(ptrows>=0, "RMatrixBDUnpackPT: PTRows<0!");
1924	if( m==0 \| n==0 \| ptrows==0 )
1925	{
1926	return;
1927	}
1928
1929	//
1930	// prepare PT
1931	//
1932	pt = new double[ptrows, n];
1933	for(i=0; i<=ptrows-1; i++)
1934	{
1935	for(j=0; j<=n-1; j++)
1936	{
1937	if( i==j )
1938	{
1939	pt[i,j] = 1;
1940	}
1941	else
1942	{
1943	pt[i,j] = 0;
1944	}
1945	}
1946	}
1947
1948	//
1949	// Calculate
1950	//
1951	rmatrixbdmultiplybyp(ref qp, m, n, ref taup, ref pt, ptrows, n, true, true);
1952	}
1953
1954
1955	/*************************************************************************
1956	Multiplication by matrix P which reduces matrix A to bidiagonal form.
1957
1958	The algorithm allows pre- or post-multiply by P or P'.
1959
1960	Input parameters:
1961	QP - matrices Q and P in compact form.
1962	Output of RMatrixBD subroutine.
1963	M - number of rows in matrix A.
1964	N - number of columns in matrix A.
1965	TAUP - scalar factors which are used to form P.
1966	Output of RMatrixBD subroutine.
1967	Z - multiplied matrix.
1968	Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
1969	ZRows - number of rows in matrix Z. If FromTheRight=False,
1970	ZRows=N, otherwise ZRows can be arbitrary.
1971	ZColumns - number of columns in matrix Z. If FromTheRight=True,
1972	ZColumns=N, otherwise ZColumns can be arbitrary.
1973	FromTheRight - pre- or post-multiply.
1974	DoTranspose - multiply by P or P'.
1975
1976	Output parameters:
1977	Z - product of Z and P.
1978	Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
1979	If ZRows=0 or ZColumns=0, the array is not modified.
1980
1981	-- ALGLIB --
1982	2005-2010
1983	Bochkanov Sergey
1984	*************************************************************************/
1985	public static void rmatrixbdmultiplybyp(ref double[,] qp,
1986	int m,
1987	int n,
1988	ref double[] taup,
1989	ref double[,] z,
1990	int zrows,
1991	int zcolumns,
1992	bool fromtheright,
1993	bool dotranspose)
1994	{
1995	int i = 0;
1996	double[] v = new double[0];
1997	double[] work = new double[0];
1998	int mx = 0;
1999	int i1 = 0;
2000	int i2 = 0;
2001	int istep = 0;
2002	int i_ = 0;
2003	int i1_ = 0;
2004
2005	if( m<=0 \| n<=0 \| zrows<=0 \| zcolumns<=0 )
2006	{
2007	return;
2008	}
2009	System.Diagnostics.Debug.Assert(fromtheright & zcolumns==n \| !fromtheright & zrows==n, "RMatrixBDMultiplyByP: incorrect Z size!");
2010
2011	//
2012	// init
2013	//
2014	mx = Math.Max(m, n);
2015	mx = Math.Max(mx, zrows);
2016	mx = Math.Max(mx, zcolumns);
2017	v = new double[mx+1];
2018	work = new double[mx+1];
2019	if( m>=n )
2020	{
2021
2022	//
2023	// setup
2024	//
2025	if( fromtheright )
2026	{
2027	i1 = n-2;
2028	i2 = 0;
2029	istep = -1;
2030	}
2031	else
2032	{
2033	i1 = 0;
2034	i2 = n-2;
2035	istep = +1;
2036	}
2037	if( !dotranspose )
2038	{
2039	i = i1;
2040	i1 = i2;
2041	i2 = i;
2042	istep = -istep;
2043	}
2044
2045	//
2046	// Process
2047	//
2048	if( n-1>0 )
2049	{
2050	i = i1;
2051	do
2052	{
2053	i1_ = (i+1) - (1);
2054	for(i_=1; i_<=n-1-i;i_++)
2055	{
2056	v[i_] = qp[i,i_+i1_];
2057	}
2058	v[1] = 1;
2059	if( fromtheright )
2060	{
2061	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i+1, n-1, ref work);
2062	}
2063	else
2064	{
2065	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i+1, n-1, 0, zcolumns-1, ref work);
2066	}
2067	i = i+istep;
2068	}
2069	while( i!=i2+istep );
2070	}
2071	}
2072	else
2073	{
2074
2075	//
2076	// setup
2077	//
2078	if( fromtheright )
2079	{
2080	i1 = m-1;
2081	i2 = 0;
2082	istep = -1;
2083	}
2084	else
2085	{
2086	i1 = 0;
2087	i2 = m-1;
2088	istep = +1;
2089	}
2090	if( !dotranspose )
2091	{
2092	i = i1;
2093	i1 = i2;
2094	i2 = i;
2095	istep = -istep;
2096	}
2097
2098	//
2099	// Process
2100	//
2101	i = i1;
2102	do
2103	{
2104	i1_ = (i) - (1);
2105	for(i_=1; i_<=n-i;i_++)
2106	{
2107	v[i_] = qp[i,i_+i1_];
2108	}
2109	v[1] = 1;
2110	if( fromtheright )
2111	{
2112	reflections.applyreflectionfromtheright(ref z, taup[i], ref v, 0, zrows-1, i, n-1, ref work);
2113	}
2114	else
2115	{
2116	reflections.applyreflectionfromtheleft(ref z, taup[i], ref v, i, n-1, 0, zcolumns-1, ref work);
2117	}
2118	i = i+istep;
2119	}
2120	while( i!=i2+istep );
2121	}
2122	}
2123
2124
2125	/*************************************************************************
2126	Unpacking of the main and secondary diagonals of bidiagonal decomposition
2127	of matrix A.
2128
2129	Input parameters:
2130	B - output of RMatrixBD subroutine.
2131	M - number of rows in matrix B.
2132	N - number of columns in matrix B.
2133
2134	Output parameters:
2135	IsUpper - True, if the matrix is upper bidiagonal.
2136	otherwise IsUpper is False.
2137	D - the main diagonal.
2138	Array whose index ranges within [0..Min(M,N)-1].
2139	E - the secondary diagonal (upper or lower, depending on
2140	the value of IsUpper).
2141	Array index ranges within [0..Min(M,N)-1], the last
2142	element is not used.
2143
2144	-- ALGLIB --
2145	2005-2010
2146	Bochkanov Sergey
2147	*************************************************************************/
2148	public static void rmatrixbdunpackdiagonals(ref double[,] b,
2149	int m,
2150	int n,
2151	ref bool isupper,
2152	ref double[] d,
2153	ref double[] e)
2154	{
2155	int i = 0;
2156
2157	isupper = m>=n;
2158	if( m<=0 \| n<=0 )
2159	{
2160	return;
2161	}
2162	if( isupper )
2163	{
2164	d = new double[n];
2165	e = new double[n];
2166	for(i=0; i<=n-2; i++)
2167	{
2168	d[i] = b[i,i];
2169	e[i] = b[i,i+1];
2170	}
2171	d[n-1] = b[n-1,n-1];
2172	}
2173	else
2174	{
2175	d = new double[m];
2176	e = new double[m];
2177	for(i=0; i<=m-2; i++)
2178	{
2179	d[i] = b[i,i];
2180	e[i] = b[i+1,i];
2181	}
2182	d[m-1] = b[m-1,m-1];
2183	}
2184	}
2185
2186
2187	/*************************************************************************
2188	Reduction of a square matrix to upper Hessenberg form: Q'AQ = H,
2189	where Q is an orthogonal matrix, H - Hessenberg matrix.
2190
2191	Input parameters:
2192	A - matrix A with elements [0..N-1, 0..N-1]
2193	N - size of matrix A.
2194
2195	Output parameters:
2196	A - matrices Q and P in compact form (see below).
2197	Tau - array of scalar factors which are used to form matrix Q.
2198	Array whose index ranges within [0..N-2]
2199
2200	Matrix H is located on the main diagonal, on the lower secondary diagonal
2201	and above the main diagonal of matrix A. The elements which are used to
2202	form matrix Q are situated in array Tau and below the lower secondary
2203	diagonal of matrix A as follows:
2204
2205	Matrix Q is represented as a product of elementary reflections
2206
2207	Q = H(0)H(2)...*H(n-2),
2208
2209	where each H(i) is given by
2210
2211	H(i) = 1 - tau * v * (v^T)
2212
2213	where tau is a scalar stored in Tau[I]; v - is a real vector,
2214	so that v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) stored in A(i+2:n-1,i).
2215
2216	-- LAPACK routine (version 3.0) --
2217	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2218	Courant Institute, Argonne National Lab, and Rice University
2219	October 31, 1992
2220	*************************************************************************/
2221	public static void rmatrixhessenberg(ref double[,] a,
2222	int n,
2223	ref double[] tau)
2224	{
2225	int i = 0;
2226	double v = 0;
2227	double[] t = new double[0];
2228	double[] work = new double[0];
2229	int i_ = 0;
2230	int i1_ = 0;
2231
2232	System.Diagnostics.Debug.Assert(n>=0, "RMatrixHessenberg: incorrect N!");
2233
2234	//
2235	// Quick return if possible
2236	//
2237	if( n<=1 )
2238	{
2239	return;
2240	}
2241	tau = new double[n-2+1];
2242	t = new double[n+1];
2243	work = new double[n-1+1];
2244	for(i=0; i<=n-2; i++)
2245	{
2246
2247	//
2248	// Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
2249	//
2250	i1_ = (i+1) - (1);
2251	for(i_=1; i_<=n-i-1;i_++)
2252	{
2253	t[i_] = a[i_+i1_,i];
2254	}
2255	reflections.generatereflection(ref t, n-i-1, ref v);
2256	i1_ = (1) - (i+1);
2257	for(i_=i+1; i_<=n-1;i_++)
2258	{
2259	a[i_,i] = t[i_+i1_];
2260	}
2261	tau[i] = v;
2262	t[1] = 1;
2263
2264	//
2265	// Apply H(i) to A(1:ihi,i+1:ihi) from the right
2266	//
2267	reflections.applyreflectionfromtheright(ref a, v, ref t, 0, n-1, i+1, n-1, ref work);
2268
2269	//
2270	// Apply H(i) to A(i+1:ihi,i+1:n) from the left
2271	//
2272	reflections.applyreflectionfromtheleft(ref a, v, ref t, i+1, n-1, i+1, n-1, ref work);
2273	}
2274	}
2275
2276
2277	/*************************************************************************
2278	Unpacking matrix Q which reduces matrix A to upper Hessenberg form
2279
2280	Input parameters:
2281	A - output of RMatrixHessenberg subroutine.
2282	N - size of matrix A.
2283	Tau - scalar factors which are used to form Q.
2284	Output of RMatrixHessenberg subroutine.
2285
2286	Output parameters:
2287	Q - matrix Q.
2288	Array whose indexes range within [0..N-1, 0..N-1].
2289
2290	-- ALGLIB --
2291	2005-2010
2292	Bochkanov Sergey
2293	*************************************************************************/
2294	public static void rmatrixhessenbergunpackq(ref double[,] a,
2295	int n,
2296	ref double[] tau,
2297	ref double[,] q)
2298	{
2299	int i = 0;
2300	int j = 0;
2301	double[] v = new double[0];
2302	double[] work = new double[0];
2303	int i_ = 0;
2304	int i1_ = 0;
2305
2306	if( n==0 )
2307	{
2308	return;
2309	}
2310
2311	//
2312	// init
2313	//
2314	q = new double[n-1+1, n-1+1];
2315	v = new double[n-1+1];
2316	work = new double[n-1+1];
2317	for(i=0; i<=n-1; i++)
2318	{
2319	for(j=0; j<=n-1; j++)
2320	{
2321	if( i==j )
2322	{
2323	q[i,j] = 1;
2324	}
2325	else
2326	{
2327	q[i,j] = 0;
2328	}
2329	}
2330	}
2331
2332	//
2333	// unpack Q
2334	//
2335	for(i=0; i<=n-2; i++)
2336	{
2337
2338	//
2339	// Apply H(i)
2340	//
2341	i1_ = (i+1) - (1);
2342	for(i_=1; i_<=n-i-1;i_++)
2343	{
2344	v[i_] = a[i_+i1_,i];
2345	}
2346	v[1] = 1;
2347	reflections.applyreflectionfromtheright(ref q, tau[i], ref v, 0, n-1, i+1, n-1, ref work);
2348	}
2349	}
2350
2351
2352	/*************************************************************************
2353	Unpacking matrix H (the result of matrix A reduction to upper Hessenberg form)
2354
2355	Input parameters:
2356	A - output of RMatrixHessenberg subroutine.
2357	N - size of matrix A.
2358
2359	Output parameters:
2360	H - matrix H. Array whose indexes range within [0..N-1, 0..N-1].
2361
2362	-- ALGLIB --
2363	2005-2010
2364	Bochkanov Sergey
2365	*************************************************************************/
2366	public static void rmatrixhessenbergunpackh(ref double[,] a,
2367	int n,
2368	ref double[,] h)
2369	{
2370	int i = 0;
2371	int j = 0;
2372	double[] v = new double[0];
2373	double[] work = new double[0];
2374	int i_ = 0;
2375
2376	if( n==0 )
2377	{
2378	return;
2379	}
2380	h = new double[n-1+1, n-1+1];
2381	for(i=0; i<=n-1; i++)
2382	{
2383	for(j=0; j<=i-2; j++)
2384	{
2385	h[i,j] = 0;
2386	}
2387	j = Math.Max(0, i-1);
2388	for(i_=j; i_<=n-1;i_++)
2389	{
2390	h[i,i_] = a[i,i_];
2391	}
2392	}
2393	}
2394
2395
2396	/*************************************************************************
2397	Reduction of a symmetric matrix which is given by its higher or lower
2398	triangular part to a tridiagonal matrix using orthogonal similarity
2399	transformation: Q'AQ=T.
2400
2401	Input parameters:
2402	A - matrix to be transformed
2403	array with elements [0..N-1, 0..N-1].
2404	N - size of matrix A.
2405	IsUpper - storage format. If IsUpper = True, then matrix A is given
2406	by its upper triangle, and the lower triangle is not used
2407	and not modified by the algorithm, and vice versa
2408	if IsUpper = False.
2409
2410	Output parameters:
2411	A - matrices T and Q in compact form (see lower)
2412	Tau - array of factors which are forming matrices H(i)
2413	array with elements [0..N-2].
2414	D - main diagonal of symmetric matrix T.
2415	array with elements [0..N-1].
2416	E - secondary diagonal of symmetric matrix T.
2417	array with elements [0..N-2].
2418
2419
2420	If IsUpper=True, the matrix Q is represented as a product of elementary
2421	reflectors
2422
2423	Q = H(n-2) . . . H(2) H(0).
2424
2425	Each H(i) has the form
2426
2427	H(i) = I - tau * v * v'
2428
2429	where tau is a real scalar, and v is a real vector with
2430	v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in
2431	A(0:i-1,i+1), and tau in TAU(i).
2432
2433	If IsUpper=False, the matrix Q is represented as a product of elementary
2434	reflectors
2435
2436	Q = H(0) H(2) . . . H(n-2).
2437
2438	Each H(i) has the form
2439
2440	H(i) = I - tau * v * v'
2441
2442	where tau is a real scalar, and v is a real vector with
2443	v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i),
2444	and tau in TAU(i).
2445
2446	The contents of A on exit are illustrated by the following examples
2447	with n = 5:
2448
2449	if UPLO = 'U': if UPLO = 'L':
2450
2451	( d e v1 v2 v3 ) ( d )
2452	( d e v2 v3 ) ( e d )
2453	( d e v3 ) ( v0 e d )
2454	( d e ) ( v0 v1 e d )
2455	( d ) ( v0 v1 v2 e d )
2456
2457	where d and e denote diagonal and off-diagonal elements of T, and vi
2458	denotes an element of the vector defining H(i).
2459
2460	-- LAPACK routine (version 3.0) --
2461	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2462	Courant Institute, Argonne National Lab, and Rice University
2463	October 31, 1992
2464	*************************************************************************/
2465	public static void smatrixtd(ref double[,] a,
2466	int n,
2467	bool isupper,
2468	ref double[] tau,
2469	ref double[] d,
2470	ref double[] e)
2471	{
2472	int i = 0;
2473	double alpha = 0;
2474	double taui = 0;
2475	double v = 0;
2476	double[] t = new double[0];
2477	double[] t2 = new double[0];
2478	double[] t3 = new double[0];
2479	int i_ = 0;
2480	int i1_ = 0;
2481
2482	if( n<=0 )
2483	{
2484	return;
2485	}
2486	t = new double[n+1];
2487	t2 = new double[n+1];
2488	t3 = new double[n+1];
2489	if( n>1 )
2490	{
2491	tau = new double[n-2+1];
2492	}
2493	d = new double[n-1+1];
2494	if( n>1 )
2495	{
2496	e = new double[n-2+1];
2497	}
2498	if( isupper )
2499	{
2500
2501	//
2502	// Reduce the upper triangle of A
2503	//
2504	for(i=n-2; i>=0; i--)
2505	{
2506
2507	//
2508	// Generate elementary reflector H() = E - tau * v * v'
2509	//
2510	if( i>=1 )
2511	{
2512	i1_ = (0) - (2);
2513	for(i_=2; i_<=i+1;i_++)
2514	{
2515	t[i_] = a[i_+i1_,i+1];
2516	}
2517	}
2518	t[1] = a[i,i+1];
2519	reflections.generatereflection(ref t, i+1, ref taui);
2520	if( i>=1 )
2521	{
2522	i1_ = (2) - (0);
2523	for(i_=0; i_<=i-1;i_++)
2524	{
2525	a[i_,i+1] = t[i_+i1_];
2526	}
2527	}
2528	a[i,i+1] = t[1];
2529	e[i] = a[i,i+1];
2530	if( (double)(taui)!=(double)(0) )
2531	{
2532
2533	//
2534	// Apply H from both sides to A
2535	//
2536	a[i,i+1] = 1;
2537
2538	//
2539	// Compute x := tau * A * v storing x in TAU
2540	//
2541	i1_ = (0) - (1);
2542	for(i_=1; i_<=i+1;i_++)
2543	{
2544	t[i_] = a[i_+i1_,i+1];
2545	}
2546	sblas.symmetricmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t3);
2547	i1_ = (1) - (0);
2548	for(i_=0; i_<=i;i_++)
2549	{
2550	tau[i_] = t3[i_+i1_];
2551	}
2552
2553	//
2554	// Compute w := x - 1/2 * tau * (x'v) v
2555	//
2556	v = 0.0;
2557	for(i_=0; i_<=i;i_++)
2558	{
2559	v += tau[i_]*a[i_,i+1];
2560	}
2561	alpha = -(0.5tauiv);
2562	for(i_=0; i_<=i;i_++)
2563	{
2564	tau[i_] = tau[i_] + alpha*a[i_,i+1];
2565	}
2566
2567	//
2568	// Apply the transformation as a rank-2 update:
2569	// A := A - v * w' - w * v'
2570	//
2571	i1_ = (0) - (1);
2572	for(i_=1; i_<=i+1;i_++)
2573	{
2574	t[i_] = a[i_+i1_,i+1];
2575	}
2576	i1_ = (0) - (1);
2577	for(i_=1; i_<=i+1;i_++)
2578	{
2579	t3[i_] = tau[i_+i1_];
2580	}
2581	sblas.symmetricrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1);
2582	a[i,i+1] = e[i];
2583	}
2584	d[i+1] = a[i+1,i+1];
2585	tau[i] = taui;
2586	}
2587	d[0] = a[0,0];
2588	}
2589	else
2590	{
2591
2592	//
2593	// Reduce the lower triangle of A
2594	//
2595	for(i=0; i<=n-2; i++)
2596	{
2597
2598	//
2599	// Generate elementary reflector H = E - tau * v * v'
2600	//
2601	i1_ = (i+1) - (1);
2602	for(i_=1; i_<=n-i-1;i_++)
2603	{
2604	t[i_] = a[i_+i1_,i];
2605	}
2606	reflections.generatereflection(ref t, n-i-1, ref taui);
2607	i1_ = (1) - (i+1);
2608	for(i_=i+1; i_<=n-1;i_++)
2609	{
2610	a[i_,i] = t[i_+i1_];
2611	}
2612	e[i] = a[i+1,i];
2613	if( (double)(taui)!=(double)(0) )
2614	{
2615
2616	//
2617	// Apply H from both sides to A
2618	//
2619	a[i+1,i] = 1;
2620
2621	//
2622	// Compute x := tau * A * v storing y in TAU
2623	//
2624	i1_ = (i+1) - (1);
2625	for(i_=1; i_<=n-i-1;i_++)
2626	{
2627	t[i_] = a[i_+i1_,i];
2628	}
2629	sblas.symmetricmatrixvectormultiply(ref a, isupper, i+1, n-1, ref t, taui, ref t2);
2630	i1_ = (1) - (i);
2631	for(i_=i; i_<=n-2;i_++)
2632	{
2633	tau[i_] = t2[i_+i1_];
2634	}
2635
2636	//
2637	// Compute w := x - 1/2 * tau * (x'v) v
2638	//
2639	i1_ = (i+1)-(i);
2640	v = 0.0;
2641	for(i_=i; i_<=n-2;i_++)
2642	{
2643	v += tau[i_]*a[i_+i1_,i];
2644	}
2645	alpha = -(0.5tauiv);
2646	i1_ = (i+1) - (i);
2647	for(i_=i; i_<=n-2;i_++)
2648	{
2649	tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
2650	}
2651
2652	//
2653	// Apply the transformation as a rank-2 update:
2654	// A := A - v * w' - w * v'
2655	//
2656	//
2657	i1_ = (i+1) - (1);
2658	for(i_=1; i_<=n-i-1;i_++)
2659	{
2660	t[i_] = a[i_+i1_,i];
2661	}
2662	i1_ = (i) - (1);
2663	for(i_=1; i_<=n-i-1;i_++)
2664	{
2665	t2[i_] = tau[i_+i1_];
2666	}
2667	sblas.symmetricrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1);
2668	a[i+1,i] = e[i];
2669	}
2670	d[i] = a[i,i];
2671	tau[i] = taui;
2672	}
2673	d[n-1] = a[n-1,n-1];
2674	}
2675	}
2676
2677
2678	/*************************************************************************
2679	Unpacking matrix Q which reduces symmetric matrix to a tridiagonal
2680	form.
2681
2682	Input parameters:
2683	A - the result of a SMatrixTD subroutine
2684	N - size of matrix A.
2685	IsUpper - storage format (a parameter of SMatrixTD subroutine)
2686	Tau - the result of a SMatrixTD subroutine
2687
2688	Output parameters:
2689	Q - transformation matrix.
2690	array with elements [0..N-1, 0..N-1].
2691
2692	-- ALGLIB --
2693	Copyright 2005-2010 by Bochkanov Sergey
2694	*************************************************************************/
2695	public static void smatrixtdunpackq(ref double[,] a,
2696	int n,
2697	bool isupper,
2698	ref double[] tau,
2699	ref double[,] q)
2700	{
2701	int i = 0;
2702	int j = 0;
2703	double[] v = new double[0];
2704	double[] work = new double[0];
2705	int i_ = 0;
2706	int i1_ = 0;
2707
2708	if( n==0 )
2709	{
2710	return;
2711	}
2712
2713	//
2714	// init
2715	//
2716	q = new double[n-1+1, n-1+1];
2717	v = new double[n+1];
2718	work = new double[n-1+1];
2719	for(i=0; i<=n-1; i++)
2720	{
2721	for(j=0; j<=n-1; j++)
2722	{
2723	if( i==j )
2724	{
2725	q[i,j] = 1;
2726	}
2727	else
2728	{
2729	q[i,j] = 0;
2730	}
2731	}
2732	}
2733
2734	//
2735	// unpack Q
2736	//
2737	if( isupper )
2738	{
2739	for(i=0; i<=n-2; i++)
2740	{
2741
2742	//
2743	// Apply H(i)
2744	//
2745	i1_ = (0) - (1);
2746	for(i_=1; i_<=i+1;i_++)
2747	{
2748	v[i_] = a[i_+i1_,i+1];
2749	}
2750	v[i+1] = 1;
2751	reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, 0, i, 0, n-1, ref work);
2752	}
2753	}
2754	else
2755	{
2756	for(i=n-2; i>=0; i--)
2757	{
2758
2759	//
2760	// Apply H(i)
2761	//
2762	i1_ = (i+1) - (1);
2763	for(i_=1; i_<=n-i-1;i_++)
2764	{
2765	v[i_] = a[i_+i1_,i];
2766	}
2767	v[1] = 1;
2768	reflections.applyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work);
2769	}
2770	}
2771	}
2772
2773
2774	/*************************************************************************
2775	Reduction of a Hermitian matrix which is given by its higher or lower
2776	triangular part to a real tridiagonal matrix using unitary similarity
2777	transformation: Q'AQ = T.
2778
2779	Input parameters:
2780	A - matrix to be transformed
2781	array with elements [0..N-1, 0..N-1].
2782	N - size of matrix A.
2783	IsUpper - storage format. If IsUpper = True, then matrix A is given
2784	by its upper triangle, and the lower triangle is not used
2785	and not modified by the algorithm, and vice versa
2786	if IsUpper = False.
2787
2788	Output parameters:
2789	A - matrices T and Q in compact form (see lower)
2790	Tau - array of factors which are forming matrices H(i)
2791	array with elements [0..N-2].
2792	D - main diagonal of real symmetric matrix T.
2793	array with elements [0..N-1].
2794	E - secondary diagonal of real symmetric matrix T.
2795	array with elements [0..N-2].
2796
2797
2798	If IsUpper=True, the matrix Q is represented as a product of elementary
2799	reflectors
2800
2801	Q = H(n-2) . . . H(2) H(0).
2802
2803	Each H(i) has the form
2804
2805	H(i) = I - tau * v * v'
2806
2807	where tau is a complex scalar, and v is a complex vector with
2808	v(i+1:n-1) = 0, v(i) = 1, v(0:i-1) is stored on exit in
2809	A(0:i-1,i+1), and tau in TAU(i).
2810
2811	If IsUpper=False, the matrix Q is represented as a product of elementary
2812	reflectors
2813
2814	Q = H(0) H(2) . . . H(n-2).
2815
2816	Each H(i) has the form
2817
2818	H(i) = I - tau * v * v'
2819
2820	where tau is a complex scalar, and v is a complex vector with
2821	v(0:i) = 0, v(i+1) = 1, v(i+2:n-1) is stored on exit in A(i+2:n-1,i),
2822	and tau in TAU(i).
2823
2824	The contents of A on exit are illustrated by the following examples
2825	with n = 5:
2826
2827	if UPLO = 'U': if UPLO = 'L':
2828
2829	( d e v1 v2 v3 ) ( d )
2830	( d e v2 v3 ) ( e d )
2831	( d e v3 ) ( v0 e d )
2832	( d e ) ( v0 v1 e d )
2833	( d ) ( v0 v1 v2 e d )
2834
2835	where d and e denote diagonal and off-diagonal elements of T, and vi
2836	denotes an element of the vector defining H(i).
2837
2838	-- LAPACK routine (version 3.0) --
2839	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2840	Courant Institute, Argonne National Lab, and Rice University
2841	October 31, 1992
2842	*************************************************************************/
2843	public static void hmatrixtd(ref AP.Complex[,] a,
2844	int n,
2845	bool isupper,
2846	ref AP.Complex[] tau,
2847	ref double[] d,
2848	ref double[] e)
2849	{
2850	int i = 0;
2851	AP.Complex alpha = 0;
2852	AP.Complex taui = 0;
2853	AP.Complex v = 0;
2854	AP.Complex[] t = new AP.Complex[0];
2855	AP.Complex[] t2 = new AP.Complex[0];
2856	AP.Complex[] t3 = new AP.Complex[0];
2857	int i_ = 0;
2858	int i1_ = 0;
2859
2860	if( n<=0 )
2861	{
2862	return;
2863	}
2864	for(i=0; i<=n-1; i++)
2865	{
2866	System.Diagnostics.Debug.Assert((double)(a[i,i].y)==(double)(0));
2867	}
2868	if( n>1 )
2869	{
2870	tau = new AP.Complex[n-2+1];
2871	e = new double[n-2+1];
2872	}
2873	d = new double[n-1+1];
2874	t = new AP.Complex[n-1+1];
2875	t2 = new AP.Complex[n-1+1];
2876	t3 = new AP.Complex[n-1+1];
2877	if( isupper )
2878	{
2879
2880	//
2881	// Reduce the upper triangle of A
2882	//
2883	a[n-1,n-1] = a[n-1,n-1].x;
2884	for(i=n-2; i>=0; i--)
2885	{
2886
2887	//
2888	// Generate elementary reflector H = I+1 - tau * v * v'
2889	//
2890	alpha = a[i,i+1];
2891	t[1] = alpha;
2892	if( i>=1 )
2893	{
2894	i1_ = (0) - (2);
2895	for(i_=2; i_<=i+1;i_++)
2896	{
2897	t[i_] = a[i_+i1_,i+1];
2898	}
2899	}
2900	creflections.complexgeneratereflection(ref t, i+1, ref taui);
2901	if( i>=1 )
2902	{
2903	i1_ = (2) - (0);
2904	for(i_=0; i_<=i-1;i_++)
2905	{
2906	a[i_,i+1] = t[i_+i1_];
2907	}
2908	}
2909	alpha = t[1];
2910	e[i] = alpha.x;
2911	if( taui!=0 )
2912	{
2913
2914	//
2915	// Apply H(I+1) from both sides to A
2916	//
2917	a[i,i+1] = 1;
2918
2919	//
2920	// Compute x := tau * A * v storing x in TAU
2921	//
2922	i1_ = (0) - (1);
2923	for(i_=1; i_<=i+1;i_++)
2924	{
2925	t[i_] = a[i_+i1_,i+1];
2926	}
2927	hblas.hermitianmatrixvectormultiply(ref a, isupper, 0, i, ref t, taui, ref t2);
2928	i1_ = (1) - (0);
2929	for(i_=0; i_<=i;i_++)
2930	{
2931	tau[i_] = t2[i_+i1_];
2932	}
2933
2934	//
2935	// Compute w := x - 1/2 * tau * (x'v) v
2936	//
2937	v = 0.0;
2938	for(i_=0; i_<=i;i_++)
2939	{
2940	v += AP.Math.Conj(tau[i_])*a[i_,i+1];
2941	}
2942	alpha = -(0.5tauiv);
2943	for(i_=0; i_<=i;i_++)
2944	{
2945	tau[i_] = tau[i_] + alpha*a[i_,i+1];
2946	}
2947
2948	//
2949	// Apply the transformation as a rank-2 update:
2950	// A := A - v * w' - w * v'
2951	//
2952	i1_ = (0) - (1);
2953	for(i_=1; i_<=i+1;i_++)
2954	{
2955	t[i_] = a[i_+i1_,i+1];
2956	}
2957	i1_ = (0) - (1);
2958	for(i_=1; i_<=i+1;i_++)
2959	{
2960	t3[i_] = tau[i_+i1_];
2961	}
2962	hblas.hermitianrank2update(ref a, isupper, 0, i, ref t, ref t3, ref t2, -1);
2963	}
2964	else
2965	{
2966	a[i,i] = a[i,i].x;
2967	}
2968	a[i,i+1] = e[i];
2969	d[i+1] = a[i+1,i+1].x;
2970	tau[i] = taui;
2971	}
2972	d[0] = a[0,0].x;
2973	}
2974	else
2975	{
2976
2977	//
2978	// Reduce the lower triangle of A
2979	//
2980	a[0,0] = a[0,0].x;
2981	for(i=0; i<=n-2; i++)
2982	{
2983
2984	//
2985	// Generate elementary reflector H = I - tau * v * v'
2986	//
2987	i1_ = (i+1) - (1);
2988	for(i_=1; i_<=n-i-1;i_++)
2989	{
2990	t[i_] = a[i_+i1_,i];
2991	}
2992	creflections.complexgeneratereflection(ref t, n-i-1, ref taui);
2993	i1_ = (1) - (i+1);
2994	for(i_=i+1; i_<=n-1;i_++)
2995	{
2996	a[i_,i] = t[i_+i1_];
2997	}
2998	e[i] = a[i+1,i].x;
2999	if( taui!=0 )
3000	{
3001
3002	//
3003	// Apply H(i) from both sides to A(i+1:n,i+1:n)
3004	//
3005	a[i+1,i] = 1;
3006
3007	//
3008	// Compute x := tau * A * v storing y in TAU
3009	//
3010	i1_ = (i+1) - (1);
3011	for(i_=1; i_<=n-i-1;i_++)
3012	{
3013	t[i_] = a[i_+i1_,i];
3014	}
3015	hblas.hermitianmatrixvectormultiply(ref a, isupper, i+1, n-1, ref t, taui, ref t2);
3016	i1_ = (1) - (i);
3017	for(i_=i; i_<=n-2;i_++)
3018	{
3019	tau[i_] = t2[i_+i1_];
3020	}
3021
3022	//
3023	// Compute w := x - 1/2 * tau * (x'v) v
3024	//
3025	i1_ = (i+1)-(i);
3026	v = 0.0;
3027	for(i_=i; i_<=n-2;i_++)
3028	{
3029	v += AP.Math.Conj(tau[i_])*a[i_+i1_,i];
3030	}
3031	alpha = -(0.5tauiv);
3032	i1_ = (i+1) - (i);
3033	for(i_=i; i_<=n-2;i_++)
3034	{
3035	tau[i_] = tau[i_] + alpha*a[i_+i1_,i];
3036	}
3037
3038	//
3039	// Apply the transformation as a rank-2 update:
3040	// A := A - v * w' - w * v'
3041	//
3042	i1_ = (i+1) - (1);
3043	for(i_=1; i_<=n-i-1;i_++)
3044	{
3045	t[i_] = a[i_+i1_,i];
3046	}
3047	i1_ = (i) - (1);
3048	for(i_=1; i_<=n-i-1;i_++)
3049	{
3050	t2[i_] = tau[i_+i1_];
3051	}
3052	hblas.hermitianrank2update(ref a, isupper, i+1, n-1, ref t, ref t2, ref t3, -1);
3053	}
3054	else
3055	{
3056	a[i+1,i+1] = a[i+1,i+1].x;
3057	}
3058	a[i+1,i] = e[i];
3059	d[i] = a[i,i].x;
3060	tau[i] = taui;
3061	}
3062	d[n-1] = a[n-1,n-1].x;
3063	}
3064	}
3065
3066
3067	/*************************************************************************
3068	Unpacking matrix Q which reduces a Hermitian matrix to a real tridiagonal
3069	form.
3070
3071	Input parameters:
3072	A - the result of a HMatrixTD subroutine
3073	N - size of matrix A.
3074	IsUpper - storage format (a parameter of HMatrixTD subroutine)
3075	Tau - the result of a HMatrixTD subroutine
3076
3077	Output parameters:
3078	Q - transformation matrix.
3079	array with elements [0..N-1, 0..N-1].
3080
3081	-- ALGLIB --
3082	Copyright 2005-2010 by Bochkanov Sergey
3083	*************************************************************************/
3084	public static void hmatrixtdunpackq(ref AP.Complex[,] a,
3085	int n,
3086	bool isupper,
3087	ref AP.Complex[] tau,
3088	ref AP.Complex[,] q)
3089	{
3090	int i = 0;
3091	int j = 0;
3092	AP.Complex[] v = new AP.Complex[0];
3093	AP.Complex[] work = new AP.Complex[0];
3094	int i_ = 0;
3095	int i1_ = 0;
3096
3097	if( n==0 )
3098	{
3099	return;
3100	}
3101
3102	//
3103	// init
3104	//
3105	q = new AP.Complex[n-1+1, n-1+1];
3106	v = new AP.Complex[n+1];
3107	work = new AP.Complex[n-1+1];
3108	for(i=0; i<=n-1; i++)
3109	{
3110	for(j=0; j<=n-1; j++)
3111	{
3112	if( i==j )
3113	{
3114	q[i,j] = 1;
3115	}
3116	else
3117	{
3118	q[i,j] = 0;
3119	}
3120	}
3121	}
3122
3123	//
3124	// unpack Q
3125	//
3126	if( isupper )
3127	{
3128	for(i=0; i<=n-2; i++)
3129	{
3130
3131	//
3132	// Apply H(i)
3133	//
3134	i1_ = (0) - (1);
3135	for(i_=1; i_<=i+1;i_++)
3136	{
3137	v[i_] = a[i_+i1_,i+1];
3138	}
3139	v[i+1] = 1;
3140	creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, 0, i, 0, n-1, ref work);
3141	}
3142	}
3143	else
3144	{
3145	for(i=n-2; i>=0; i--)
3146	{
3147
3148	//
3149	// Apply H(i)
3150	//
3151	i1_ = (i+1) - (1);
3152	for(i_=1; i_<=n-i-1;i_++)
3153	{
3154	v[i_] = a[i_+i1_,i];
3155	}
3156	v[1] = 1;
3157	creflections.complexapplyreflectionfromtheleft(ref q, tau[i], ref v, i+1, n-1, 0, n-1, ref work);
3158	}
3159	}
3160	}
3161
3162
3163	/*************************************************************************
3164	Base case for real QR
3165
3166	-- LAPACK routine (version 3.0) --
3167	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3168	Courant Institute, Argonne National Lab, and Rice University
3169	September 30, 1994.
3170	Sergey Bochkanov, ALGLIB project, translation from FORTRAN to
3171	pseudocode, 2007-2010.
3172	*************************************************************************/
3173	private static void rmatrixqrbasecase(ref double[,] a,
3174	int m,
3175	int n,
3176	ref double[] work,
3177	ref double[] t,
3178	ref double[] tau)
3179	{
3180	int i = 0;
3181	int k = 0;
3182	int minmn = 0;
3183	double tmp = 0;
3184	int i_ = 0;
3185	int i1_ = 0;
3186
3187	minmn = Math.Min(m, n);
3188
3189	//
3190	// Test the input arguments
3191	//
3192	k = minmn;
3193	for(i=0; i<=k-1; i++)
3194	{
3195
3196	//
3197	// Generate elementary reflector H(i) to annihilate A(i+1:m,i)
3198	//
3199	i1_ = (i) - (1);
3200	for(i_=1; i_<=m-i;i_++)
3201	{
3202	t[i_] = a[i_+i1_,i];
3203	}
3204	reflections.generatereflection(ref t, m-i, ref tmp);
3205	tau[i] = tmp;
3206	i1_ = (1) - (i);
3207	for(i_=i; i_<=m-1;i_++)
3208	{
3209	a[i_,i] = t[i_+i1_];
3210	}
3211	t[1] = 1;
3212	if( i<n )
3213	{
3214
3215	//
3216	// Apply H(i) to A(i:m-1,i+1:n-1) from the left
3217	//
3218	reflections.applyreflectionfromtheleft(ref a, tau[i], ref t, i, m-1, i+1, n-1, ref work);
3219	}
3220	}
3221	}
3222
3223
3224	/*************************************************************************
3225	Base case for real LQ
3226
3227	-- LAPACK routine (version 3.0) --
3228	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3229	Courant Institute, Argonne National Lab, and Rice University
3230	September 30, 1994.
3231	Sergey Bochkanov, ALGLIB project, translation from FORTRAN to
3232	pseudocode, 2007-2010.
3233	*************************************************************************/
3234	private static void rmatrixlqbasecase(ref double[,] a,
3235	int m,
3236	int n,
3237	ref double[] work,
3238	ref double[] t,
3239	ref double[] tau)
3240	{
3241	int i = 0;
3242	int k = 0;
3243	int minmn = 0;
3244	double tmp = 0;
3245	int i_ = 0;
3246	int i1_ = 0;
3247
3248	minmn = Math.Min(m, n);
3249	k = Math.Min(m, n);
3250	for(i=0; i<=k-1; i++)
3251	{
3252
3253	//
3254	// Generate elementary reflector H(i) to annihilate A(i,i+1:n-1)
3255	//
3256	i1_ = (i) - (1);
3257	for(i_=1; i_<=n-i;i_++)
3258	{
3259	t[i_] = a[i,i_+i1_];
3260	}
3261	reflections.generatereflection(ref t, n-i, ref tmp);
3262	tau[i] = tmp;
3263	i1_ = (1) - (i);
3264	for(i_=i; i_<=n-1;i_++)
3265	{
3266	a[i,i_] = t[i_+i1_];
3267	}
3268	t[1] = 1;
3269	if( i<n )
3270	{
3271
3272	//
3273	// Apply H(i) to A(i+1:m,i:n) from the right
3274	//
3275	reflections.applyreflectionfromtheright(ref a, tau[i], ref t, i+1, m-1, i, n-1, ref work);
3276	}
3277	}
3278	}
3279
3280
3281	/*************************************************************************
3282	Base case for complex QR
3283
3284	-- LAPACK routine (version 3.0) --
3285	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3286	Courant Institute, Argonne National Lab, and Rice University
3287	September 30, 1994.
3288	Sergey Bochkanov, ALGLIB project, translation from FORTRAN to
3289	pseudocode, 2007-2010.
3290	*************************************************************************/
3291	private static void cmatrixqrbasecase(ref AP.Complex[,] a,
3292	int m,
3293	int n,
3294	ref AP.Complex[] work,
3295	ref AP.Complex[] t,
3296	ref AP.Complex[] tau)
3297	{
3298	int i = 0;
3299	int k = 0;
3300	int mmi = 0;
3301	int minmn = 0;
3302	AP.Complex tmp = 0;
3303	int i_ = 0;
3304	int i1_ = 0;
3305
3306	minmn = Math.Min(m, n);
3307	if( minmn<=0 )
3308	{
3309	return;
3310	}
3311
3312	//
3313	// Test the input arguments
3314	//
3315	k = Math.Min(m, n);
3316	for(i=0; i<=k-1; i++)
3317	{
3318
3319	//
3320	// Generate elementary reflector H(i) to annihilate A(i+1:m,i)
3321	//
3322	mmi = m-i;
3323	i1_ = (i) - (1);
3324	for(i_=1; i_<=mmi;i_++)
3325	{
3326	t[i_] = a[i_+i1_,i];
3327	}
3328	creflections.complexgeneratereflection(ref t, mmi, ref tmp);
3329	tau[i] = tmp;
3330	i1_ = (1) - (i);
3331	for(i_=i; i_<=m-1;i_++)
3332	{
3333	a[i_,i] = t[i_+i1_];
3334	}
3335	t[1] = 1;
3336	if( i<n-1 )
3337	{
3338
3339	//
3340	// Apply H'(i) to A(i:m,i+1:n) from the left
3341	//
3342	creflections.complexapplyreflectionfromtheleft(ref a, AP.Math.Conj(tau[i]), ref t, i, m-1, i+1, n-1, ref work);
3343	}
3344	}
3345	}
3346
3347
3348	/*************************************************************************
3349	Base case for complex LQ
3350
3351	-- LAPACK routine (version 3.0) --
3352	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
3353	Courant Institute, Argonne National Lab, and Rice University
3354	September 30, 1994.
3355	Sergey Bochkanov, ALGLIB project, translation from FORTRAN to
3356	pseudocode, 2007-2010.
3357	*************************************************************************/
3358	private static void cmatrixlqbasecase(ref AP.Complex[,] a,
3359	int m,
3360	int n,
3361	ref AP.Complex[] work,
3362	ref AP.Complex[] t,
3363	ref AP.Complex[] tau)
3364	{
3365	int i = 0;
3366	int minmn = 0;
3367	AP.Complex tmp = 0;
3368	int i_ = 0;
3369	int i1_ = 0;
3370
3371	minmn = Math.Min(m, n);
3372	if( minmn<=0 )
3373	{
3374	return;
3375	}
3376
3377	//
3378	// Test the input arguments
3379	//
3380	for(i=0; i<=minmn-1; i++)
3381	{
3382
3383	//
3384	// Generate elementary reflector H(i)
3385	//
3386	// NOTE: ComplexGenerateReflection() generates left reflector,
3387	// i.e. H which reduces x by applyiong from the left, but we
3388	// need RIGHT reflector. So we replace H=E-tauvv' by H^H,
3389	// which changes v to conj(v).
3390	//
3391	i1_ = (i) - (1);
3392	for(i_=1; i_<=n-i;i_++)
3393	{
3394	t[i_] = AP.Math.Conj(a[i,i_+i1_]);
3395	}
3396	creflections.complexgeneratereflection(ref t, n-i, ref tmp);
3397	tau[i] = tmp;
3398	i1_ = (1) - (i);
3399	for(i_=i; i_<=n-1;i_++)
3400	{
3401	a[i,i_] = AP.Math.Conj(t[i_+i1_]);
3402	}
3403	t[1] = 1;
3404	if( i<m-1 )
3405	{
3406
3407	//
3408	// Apply H'(i)
3409	//
3410	creflections.complexapplyreflectionfromtheright(ref a, tau[i], ref t, i+1, m-1, i, n-1, ref work);
3411	}
3412	}
3413	}
3414
3415
3416	/*************************************************************************
3417	Generate block reflector:
3418	* fill unused parts of reflectors matrix by zeros
3419	* fill diagonal of reflectors matrix by ones
3420	* generate triangular factor T
3421
3422	PARAMETERS:
3423	A - either LengthA*BlockSize (if ColumnwiseA) or
3424	BlockSize*LengthA (if not ColumnwiseA) matrix of
3425	elementary reflectors.
3426	Modified on exit.
3427	Tau - scalar factors
3428	ColumnwiseA - reflectors are stored in rows or in columns
3429	LengthA - length of largest reflector
3430	BlockSize - number of reflectors
3431	T - array[BlockSize,2BlockSize]. Left BlockSizeBlockSize
3432	submatrix stores triangular factor on exit.
3433	WORK - array[BlockSize]
3434
3435	-- ALGLIB routine --
3436	17.02.2010
3437	Bochkanov Sergey
3438	*************************************************************************/
3439	private static void rmatrixblockreflector(ref double[,] a,
3440	ref double[] tau,
3441	bool columnwisea,
3442	int lengtha,
3443	int blocksize,
3444	ref double[,] t,
3445	ref double[] work)
3446	{
3447	int i = 0;
3448	int j = 0;
3449	int k = 0;
3450	double v = 0;
3451	int i_ = 0;
3452	int i1_ = 0;
3453
3454
3455	//
3456	// fill beginning of new column with zeros,
3457	// load 1.0 in the first non-zero element
3458	//
3459	for(k=0; k<=blocksize-1; k++)
3460	{
3461	if( columnwisea )
3462	{
3463	for(i=0; i<=k-1; i++)
3464	{
3465	a[i,k] = 0;
3466	}
3467	}
3468	else
3469	{
3470	for(i=0; i<=k-1; i++)
3471	{
3472	a[k,i] = 0;
3473	}
3474	}
3475	a[k,k] = 1;
3476	}
3477
3478	//
3479	// Calculate Gram matrix of A
3480	//
3481	for(i=0; i<=blocksize-1; i++)
3482	{
3483	for(j=0; j<=blocksize-1; j++)
3484	{
3485	t[i,blocksize+j] = 0;
3486	}
3487	}
3488	for(k=0; k<=lengtha-1; k++)
3489	{
3490	for(j=1; j<=blocksize-1; j++)
3491	{
3492	if( columnwisea )
3493	{
3494	v = a[k,j];
3495	if( (double)(v)!=(double)(0) )
3496	{
3497	i1_ = (0) - (blocksize);
3498	for(i_=blocksize; i_<=blocksize+j-1;i_++)
3499	{
3500	t[j,i_] = t[j,i_] + v*a[k,i_+i1_];
3501	}
3502	}
3503	}
3504	else
3505	{
3506	v = a[j,k];
3507	if( (double)(v)!=(double)(0) )
3508	{
3509	i1_ = (0) - (blocksize);
3510	for(i_=blocksize; i_<=blocksize+j-1;i_++)
3511	{
3512	t[j,i_] = t[j,i_] + v*a[i_+i1_,k];
3513	}
3514	}
3515	}
3516	}
3517	}
3518
3519	//
3520	// Prepare Y (stored in TmpA) and T (stored in TmpT)
3521	//
3522	for(k=0; k<=blocksize-1; k++)
3523	{
3524
3525	//
3526	// fill non-zero part of T, use pre-calculated Gram matrix
3527	//
3528	i1_ = (blocksize) - (0);
3529	for(i_=0; i_<=k-1;i_++)
3530	{
3531	work[i_] = t[k,i_+i1_];
3532	}
3533	for(i=0; i<=k-1; i++)
3534	{
3535	v = 0.0;
3536	for(i_=i; i_<=k-1;i_++)
3537	{
3538	v += t[i,i_]*work[i_];
3539	}
3540	t[i,k] = -(tau[k]*v);
3541	}
3542	t[k,k] = -tau[k];
3543
3544	//
3545	// Rest of T is filled by zeros
3546	//
3547	for(i=k+1; i<=blocksize-1; i++)
3548	{
3549	t[i,k] = 0;
3550	}
3551	}
3552	}
3553
3554
3555	/*************************************************************************
3556	Generate block reflector (complex):
3557	* fill unused parts of reflectors matrix by zeros
3558	* fill diagonal of reflectors matrix by ones
3559	* generate triangular factor T
3560
3561
3562	-- ALGLIB routine --
3563	17.02.2010
3564	Bochkanov Sergey
3565	*************************************************************************/
3566	private static void cmatrixblockreflector(ref AP.Complex[,] a,
3567	ref AP.Complex[] tau,
3568	bool columnwisea,
3569	int lengtha,
3570	int blocksize,
3571	ref AP.Complex[,] t,
3572	ref AP.Complex[] work)
3573	{
3574	int i = 0;
3575	int k = 0;
3576	AP.Complex v = 0;
3577	int i_ = 0;
3578
3579
3580	//
3581	// Prepare Y (stored in TmpA) and T (stored in TmpT)
3582	//
3583	for(k=0; k<=blocksize-1; k++)
3584	{
3585
3586	//
3587	// fill beginning of new column with zeros,
3588	// load 1.0 in the first non-zero element
3589	//
3590	if( columnwisea )
3591	{
3592	for(i=0; i<=k-1; i++)
3593	{
3594	a[i,k] = 0;
3595	}
3596	}
3597	else
3598	{
3599	for(i=0; i<=k-1; i++)
3600	{
3601	a[k,i] = 0;
3602	}
3603	}
3604	a[k,k] = 1;
3605
3606	//
3607	// fill non-zero part of T,
3608	//
3609	for(i=0; i<=k-1; i++)
3610	{
3611	if( columnwisea )
3612	{
3613	v = 0.0;
3614	for(i_=k; i_<=lengtha-1;i_++)
3615	{
3616	v += AP.Math.Conj(a[i_,i])*a[i_,k];
3617	}
3618	}
3619	else
3620	{
3621	v = 0.0;
3622	for(i_=k; i_<=lengtha-1;i_++)
3623	{
3624	v += a[i,i_]*AP.Math.Conj(a[k,i_]);
3625	}
3626	}
3627	work[i] = v;
3628	}
3629	for(i=0; i<=k-1; i++)
3630	{
3631	v = 0.0;
3632	for(i_=i; i_<=k-1;i_++)
3633	{
3634	v += t[i,i_]*work[i_];
3635	}
3636	t[i,k] = -(tau[k]*v);
3637	}
3638	t[k,k] = -tau[k];
3639
3640	//
3641	// Rest of T is filled by zeros
3642	//
3643	for(i=k+1; i<=blocksize-1; i++)
3644	{
3645	t[i,k] = 0;
3646	}
3647	}
3648	}
3649	}
3650	}

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

Download in other formats:

Update cookies preferences