Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/matinv.cs @ 3839

Visit:

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

Line
1	/*************************************************************************
2	Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
3
4	Contributors:
5	* Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6	pseudocode.
7
8	See subroutines comments for additional copyrights.
9
10	>>> SOURCE LICENSE >>>
11	This program is free software; you can redistribute it and/or modify
12	it under the terms of the GNU General Public License as published by
13	the Free Software Foundation (www.fsf.org); either version 2 of the
14	License, or (at your option) any later version.
15
16	This program is distributed in the hope that it will be useful,
17	but WITHOUT ANY WARRANTY; without even the implied warranty of
18	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19	GNU General Public License for more details.
20
21	A copy of the GNU General Public License is available at
22	http://www.fsf.org/licensing/licenses
23
24	>>> END OF LICENSE >>>
25	*************************************************************************/
26
27	using System;
28
29	namespace alglib
30	{
31	public class matinv
32	{
33	public struct matinvreport
34	{
35	public double r1;
36	public double rinf;
37	};
38
39
40
41
42	/*************************************************************************
43	Inversion of a matrix given by its LU decomposition.
44
45	INPUT PARAMETERS:
46	A - LU decomposition of the matrix (output of RMatrixLU subroutine).
47	Pivots - table of permutations which were made during the LU decomposition
48	(the output of RMatrixLU subroutine).
49	N - size of matrix A.
50
51	OUTPUT PARAMETERS:
52	Info - return code:
53	* -3 A is singular, or VERY close to singular.
54	it is filled by zeros in such cases.
55	* -1 N<=0 was passed, or incorrect Pivots was passed
56	* 1 task is solved (but matrix A may be ill-conditioned,
57	check R1/RInf parameters for condition numbers).
58	Rep - solver report, see below for more info
59	A - inverse of matrix A.
60	Array whose indexes range within [0..N-1, 0..N-1].
61
62	SOLVER REPORT
63
64	Subroutine sets following fields of the Rep structure:
65	* R1 reciprocal of condition number: 1/cond(A), 1-norm.
66	* RInf reciprocal of condition number: 1/cond(A), inf-norm.
67
68	-- ALGLIB routine --
69	05.02.2010
70	Bochkanov Sergey
71	*************************************************************************/
72	public static void rmatrixluinverse(ref double[,] a,
73	ref int[] pivots,
74	int n,
75	ref int info,
76	ref matinvreport rep)
77	{
78	double[] work = new double[0];
79	int i = 0;
80	int j = 0;
81	int k = 0;
82	double v = 0;
83
84	info = 1;
85
86	//
87	// Quick return if possible
88	//
89	if( n==0 )
90	{
91	info = -1;
92	return;
93	}
94	for(i=0; i<=n-1; i++)
95	{
96	if( pivots[i]>n-1 \| pivots[i]<i )
97	{
98	info = -1;
99	return;
100	}
101	}
102
103	//
104	// calculate condition numbers
105	//
106	rep.r1 = rcond.rmatrixlurcond1(ref a, n);
107	rep.rinf = rcond.rmatrixlurcondinf(ref a, n);
108	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
109	{
110	for(i=0; i<=n-1; i++)
111	{
112	for(j=0; j<=n-1; j++)
113	{
114	a[i,j] = 0;
115	}
116	}
117	rep.r1 = 0;
118	rep.rinf = 0;
119	info = -3;
120	return;
121	}
122
123	//
124	// Call cache-oblivious code
125	//
126	work = new double[n];
127	rmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);
128
129	//
130	// apply permutations
131	//
132	for(i=0; i<=n-1; i++)
133	{
134	for(j=n-2; j>=0; j--)
135	{
136	k = pivots[j];
137	v = a[i,j];
138	a[i,j] = a[i,k];
139	a[i,k] = v;
140	}
141	}
142	}
143
144
145	/*************************************************************************
146	Inversion of a general matrix.
147
148	Input parameters:
149	A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
150	N - size of matrix A.
151
152	Output parameters:
153	Info - return code, same as in RMatrixLUInverse
154	Rep - solver report, same as in RMatrixLUInverse
155	A - inverse of matrix A, same as in RMatrixLUInverse
156
157	Result:
158	True, if the matrix is not singular.
159	False, if the matrix is singular.
160
161	-- ALGLIB --
162	Copyright 2005 by Bochkanov Sergey
163	*************************************************************************/
164	public static void rmatrixinverse(ref double[,] a,
165	int n,
166	ref int info,
167	ref matinvreport rep)
168	{
169	int[] pivots = new int[0];
170
171	trfac.rmatrixlu(ref a, n, n, ref pivots);
172	rmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
173	}
174
175
176	/*************************************************************************
177	Inversion of a matrix given by its LU decomposition.
178
179	INPUT PARAMETERS:
180	A - LU decomposition of the matrix (output of CMatrixLU subroutine).
181	Pivots - table of permutations which were made during the LU decomposition
182	(the output of CMatrixLU subroutine).
183	N - size of matrix A.
184
185	OUTPUT PARAMETERS:
186	Info - return code, same as in RMatrixLUInverse
187	Rep - solver report, same as in RMatrixLUInverse
188	A - inverse of matrix A, same as in RMatrixLUInverse
189
190	-- ALGLIB routine --
191	05.02.2010
192	Bochkanov Sergey
193	*************************************************************************/
194	public static void cmatrixluinverse(ref AP.Complex[,] a,
195	ref int[] pivots,
196	int n,
197	ref int info,
198	ref matinvreport rep)
199	{
200	AP.Complex[] work = new AP.Complex[0];
201	int i = 0;
202	int j = 0;
203	int k = 0;
204	AP.Complex v = 0;
205
206	info = 1;
207
208	//
209	// Quick return if possible
210	//
211	if( n==0 )
212	{
213	info = -1;
214	return;
215	}
216	for(i=0; i<=n-1; i++)
217	{
218	if( pivots[i]>n-1 \| pivots[i]<i )
219	{
220	info = -1;
221	return;
222	}
223	}
224
225	//
226	// calculate condition numbers
227	//
228	rep.r1 = rcond.cmatrixlurcond1(ref a, n);
229	rep.rinf = rcond.cmatrixlurcondinf(ref a, n);
230	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
231	{
232	for(i=0; i<=n-1; i++)
233	{
234	for(j=0; j<=n-1; j++)
235	{
236	a[i,j] = 0;
237	}
238	}
239	rep.r1 = 0;
240	rep.rinf = 0;
241	info = -3;
242	return;
243	}
244
245	//
246	// Call cache-oblivious code
247	//
248	work = new AP.Complex[n];
249	cmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);
250
251	//
252	// apply permutations
253	//
254	for(i=0; i<=n-1; i++)
255	{
256	for(j=n-2; j>=0; j--)
257	{
258	k = pivots[j];
259	v = a[i,j];
260	a[i,j] = a[i,k];
261	a[i,k] = v;
262	}
263	}
264	}
265
266
267	/*************************************************************************
268	Inversion of a general matrix.
269
270	Input parameters:
271	A - matrix, array[0..N-1,0..N-1].
272	N - size of A.
273
274	Output parameters:
275	Info - return code, same as in RMatrixLUInverse
276	Rep - solver report, same as in RMatrixLUInverse
277	A - inverse of matrix A, same as in RMatrixLUInverse
278
279	-- ALGLIB --
280	Copyright 2005 by Bochkanov Sergey
281	*************************************************************************/
282	public static void cmatrixinverse(ref AP.Complex[,] a,
283	int n,
284	ref int info,
285	ref matinvreport rep)
286	{
287	int[] pivots = new int[0];
288
289	trfac.cmatrixlu(ref a, n, n, ref pivots);
290	cmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
291	}
292
293
294	/*************************************************************************
295	Inversion of a symmetric positive definite matrix which is given
296	by Cholesky decomposition.
297
298	Input parameters:
299	A - Cholesky decomposition of the matrix to be inverted:
300	A=UU or A = LL'.
301	Output of SPDMatrixCholesky subroutine.
302	N - size of matrix A.
303	IsUpper storage format.
304	If IsUpper = True, then matrix A is given as A = U'*U
305	(matrix contains upper triangle).
306	Similarly, if IsUpper = False, then A = L*L'.
307
308	Output parameters:
309	Info - return code, same as in RMatrixLUInverse
310	Rep - solver report, same as in RMatrixLUInverse
311	A - inverse of matrix A, same as in RMatrixLUInverse
312
313	-- ALGLIB routine --
314	10.02.2010
315	Bochkanov Sergey
316	*************************************************************************/
317	public static void spdmatrixcholeskyinverse(ref double[,] a,
318	int n,
319	bool isupper,
320	ref int info,
321	ref matinvreport rep)
322	{
323	int i = 0;
324	int j = 0;
325	int k = 0;
326	double v = 0;
327	double ajj = 0;
328	double aii = 0;
329	double[] tmp = new double[0];
330	int info2 = 0;
331	matinvreport rep2 = new matinvreport();
332
333	if( n<1 )
334	{
335	info = -1;
336	return;
337	}
338	info = 1;
339
340	//
341	// calculate condition numbers
342	//
343	rep.r1 = rcond.spdmatrixcholeskyrcond(ref a, n, isupper);
344	rep.rinf = rep.r1;
345	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
346	{
347	if( isupper )
348	{
349	for(i=0; i<=n-1; i++)
350	{
351	for(j=i; j<=n-1; j++)
352	{
353	a[i,j] = 0;
354	}
355	}
356	}
357	else
358	{
359	for(i=0; i<=n-1; i++)
360	{
361	for(j=0; j<=i; j++)
362	{
363	a[i,j] = 0;
364	}
365	}
366	}
367	rep.r1 = 0;
368	rep.rinf = 0;
369	info = -3;
370	return;
371	}
372
373	//
374	// Inverse
375	//
376	tmp = new double[n];
377	spdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
378	}
379
380
381	/*************************************************************************
382	Inversion of a symmetric positive definite matrix.
383
384	Given an upper or lower triangle of a symmetric positive definite matrix,
385	the algorithm generates matrix A^-1 and saves the upper or lower triangle
386	depending on the input.
387
388	Input parameters:
389	A - matrix to be inverted (upper or lower triangle).
390	Array with elements [0..N-1,0..N-1].
391	N - size of matrix A.
392	IsUpper - storage format.
393	If IsUpper = True, then the upper triangle of matrix A is
394	given, otherwise the lower triangle is given.
395
396	Output parameters:
397	Info - return code, same as in RMatrixLUInverse
398	Rep - solver report, same as in RMatrixLUInverse
399	A - inverse of matrix A, same as in RMatrixLUInverse
400
401	-- ALGLIB routine --
402	10.02.2010
403	Bochkanov Sergey
404	*************************************************************************/
405	public static void spdmatrixinverse(ref double[,] a,
406	int n,
407	bool isupper,
408	ref int info,
409	ref matinvreport rep)
410	{
411	if( n<1 )
412	{
413	info = -1;
414	return;
415	}
416	info = 1;
417	if( trfac.spdmatrixcholesky(ref a, n, isupper) )
418	{
419	spdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
420	}
421	else
422	{
423	info = -3;
424	}
425	}
426
427
428	/*************************************************************************
429	Inversion of a Hermitian positive definite matrix which is given
430	by Cholesky decomposition.
431
432	Input parameters:
433	A - Cholesky decomposition of the matrix to be inverted:
434	A=UU or A = LL'.
435	Output of HPDMatrixCholesky subroutine.
436	N - size of matrix A.
437	IsUpper storage format.
438	If IsUpper = True, then matrix A is given as A = U'*U
439	(matrix contains upper triangle).
440	Similarly, if IsUpper = False, then A = L*L'.
441
442	Output parameters:
443	Info - return code, same as in RMatrixLUInverse
444	Rep - solver report, same as in RMatrixLUInverse
445	A - inverse of matrix A, same as in RMatrixLUInverse
446
447	-- ALGLIB routine --
448	10.02.2010
449	Bochkanov Sergey
450	*************************************************************************/
451	public static void hpdmatrixcholeskyinverse(ref AP.Complex[,] a,
452	int n,
453	bool isupper,
454	ref int info,
455	ref matinvreport rep)
456	{
457	int i = 0;
458	int j = 0;
459	int info2 = 0;
460	matinvreport rep2 = new matinvreport();
461	AP.Complex[] tmp = new AP.Complex[0];
462	AP.Complex v = 0;
463
464	if( n<1 )
465	{
466	info = -1;
467	return;
468	}
469	info = 1;
470
471	//
472	// calculate condition numbers
473	//
474	rep.r1 = rcond.hpdmatrixcholeskyrcond(ref a, n, isupper);
475	rep.rinf = rep.r1;
476	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
477	{
478	if( isupper )
479	{
480	for(i=0; i<=n-1; i++)
481	{
482	for(j=i; j<=n-1; j++)
483	{
484	a[i,j] = 0;
485	}
486	}
487	}
488	else
489	{
490	for(i=0; i<=n-1; i++)
491	{
492	for(j=0; j<=i; j++)
493	{
494	a[i,j] = 0;
495	}
496	}
497	}
498	rep.r1 = 0;
499	rep.rinf = 0;
500	info = -3;
501	return;
502	}
503
504	//
505	// Inverse
506	//
507	tmp = new AP.Complex[n];
508	hpdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
509	}
510
511
512	/*************************************************************************
513	Inversion of a Hermitian positive definite matrix.
514
515	Given an upper or lower triangle of a Hermitian positive definite matrix,
516	the algorithm generates matrix A^-1 and saves the upper or lower triangle
517	depending on the input.
518
519	Input parameters:
520	A - matrix to be inverted (upper or lower triangle).
521	Array with elements [0..N-1,0..N-1].
522	N - size of matrix A.
523	IsUpper - storage format.
524	If IsUpper = True, then the upper triangle of matrix A is
525	given, otherwise the lower triangle is given.
526
527	Output parameters:
528	Info - return code, same as in RMatrixLUInverse
529	Rep - solver report, same as in RMatrixLUInverse
530	A - inverse of matrix A, same as in RMatrixLUInverse
531
532	-- ALGLIB routine --
533	10.02.2010
534	Bochkanov Sergey
535	*************************************************************************/
536	public static void hpdmatrixinverse(ref AP.Complex[,] a,
537	int n,
538	bool isupper,
539	ref int info,
540	ref matinvreport rep)
541	{
542	if( n<1 )
543	{
544	info = -1;
545	return;
546	}
547	info = 1;
548	if( trfac.hpdmatrixcholesky(ref a, n, isupper) )
549	{
550	hpdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
551	}
552	else
553	{
554	info = -3;
555	}
556	}
557
558
559	/*************************************************************************
560	Triangular matrix inverse (real)
561
562	The subroutine inverts the following types of matrices:
563	* upper triangular
564	* upper triangular with unit diagonal
565	* lower triangular
566	* lower triangular with unit diagonal
567
568	In case of an upper (lower) triangular matrix, the inverse matrix will
569	also be upper (lower) triangular, and after the end of the algorithm, the
570	inverse matrix replaces the source matrix. The elements below (above) the
571	main diagonal are not changed by the algorithm.
572
573	If the matrix has a unit diagonal, the inverse matrix also has a unit
574	diagonal, and the diagonal elements are not passed to the algorithm.
575
576	Input parameters:
577	A - matrix, array[0..N-1, 0..N-1].
578	N - size of A.
579	IsUpper - True, if the matrix is upper triangular.
580	IsUnit - True, if the matrix has a unit diagonal.
581
582	Output parameters:
583	Info - same as for RMatrixLUInverse
584	Rep - same as for RMatrixLUInverse
585	A - same as for RMatrixLUInverse.
586
587	-- ALGLIB --
588	Copyright 05.02.2010 by Bochkanov Sergey
589	*************************************************************************/
590	public static void rmatrixtrinverse(ref double[,] a,
591	int n,
592	bool isupper,
593	bool isunit,
594	ref int info,
595	ref matinvreport rep)
596	{
597	int i = 0;
598	int j = 0;
599	double[] tmp = new double[0];
600
601	if( n<1 )
602	{
603	info = -1;
604	return;
605	}
606	info = 1;
607
608	//
609	// calculate condition numbers
610	//
611	rep.r1 = rcond.rmatrixtrrcond1(ref a, n, isupper, isunit);
612	rep.rinf = rcond.rmatrixtrrcondinf(ref a, n, isupper, isunit);
613	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
614	{
615	for(i=0; i<=n-1; i++)
616	{
617	for(j=0; j<=n-1; j++)
618	{
619	a[i,j] = 0;
620	}
621	}
622	rep.r1 = 0;
623	rep.rinf = 0;
624	info = -3;
625	return;
626	}
627
628	//
629	// Invert
630	//
631	tmp = new double[n];
632	rmatrixtrinverserec(ref a, 0, n, isupper, isunit, ref tmp, ref info, ref rep);
633	}
634
635
636	/*************************************************************************
637	Triangular matrix inverse (complex)
638
639	The subroutine inverts the following types of matrices:
640	* upper triangular
641	* upper triangular with unit diagonal
642	* lower triangular
643	* lower triangular with unit diagonal
644
645	In case of an upper (lower) triangular matrix, the inverse matrix will
646	also be upper (lower) triangular, and after the end of the algorithm, the
647	inverse matrix replaces the source matrix. The elements below (above) the
648	main diagonal are not changed by the algorithm.
649
650	If the matrix has a unit diagonal, the inverse matrix also has a unit
651	diagonal, and the diagonal elements are not passed to the algorithm.
652
653	Input parameters:
654	A - matrix, array[0..N-1, 0..N-1].
655	N - size of A.
656	IsUpper - True, if the matrix is upper triangular.
657	IsUnit - True, if the matrix has a unit diagonal.
658
659	Output parameters:
660	Info - same as for RMatrixLUInverse
661	Rep - same as for RMatrixLUInverse
662	A - same as for RMatrixLUInverse.
663
664	-- ALGLIB --
665	Copyright 05.02.2010 by Bochkanov Sergey
666	*************************************************************************/
667	public static void cmatrixtrinverse(ref AP.Complex[,] a,
668	int n,
669	bool isupper,
670	bool isunit,
671	ref int info,
672	ref matinvreport rep)
673	{
674	int i = 0;
675	int j = 0;
676	AP.Complex[] tmp = new AP.Complex[0];
677
678	if( n<1 )
679	{
680	info = -1;
681	return;
682	}
683	info = 1;
684
685	//
686	// calculate condition numbers
687	//
688	rep.r1 = rcond.cmatrixtrrcond1(ref a, n, isupper, isunit);
689	rep.rinf = rcond.cmatrixtrrcondinf(ref a, n, isupper, isunit);
690	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
691	{
692	for(i=0; i<=n-1; i++)
693	{
694	for(j=0; j<=n-1; j++)
695	{
696	a[i,j] = 0;
697	}
698	}
699	rep.r1 = 0;
700	rep.rinf = 0;
701	info = -3;
702	return;
703	}
704
705	//
706	// Invert
707	//
708	tmp = new AP.Complex[n];
709	cmatrixtrinverserec(ref a, 0, n, isupper, isunit, ref tmp, ref info, ref rep);
710	}
711
712
713	/*************************************************************************
714	Triangular matrix inversion, recursive subroutine
715
716	-- ALGLIB --
717	05.02.2010, Bochkanov Sergey.
718	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
719	Courant Institute, Argonne National Lab, and Rice University
720	February 29, 1992.
721	*************************************************************************/
722	private static void rmatrixtrinverserec(ref double[,] a,
723	int offs,
724	int n,
725	bool isupper,
726	bool isunit,
727	ref double[] tmp,
728	ref int info,
729	ref matinvreport rep)
730	{
731	int n1 = 0;
732	int n2 = 0;
733	int i = 0;
734	int j = 0;
735	double v = 0;
736	double ajj = 0;
737	int i_ = 0;
738	int i1_ = 0;
739
740	if( n<1 )
741	{
742	info = -1;
743	return;
744	}
745
746	//
747	// Base case
748	//
749	if( n<=ablas.ablasblocksize(ref a) )
750	{
751	if( isupper )
752	{
753
754	//
755	// Compute inverse of upper triangular matrix.
756	//
757	for(j=0; j<=n-1; j++)
758	{
759	if( !isunit )
760	{
761	if( (double)(a[offs+j,offs+j])==(double)(0) )
762	{
763	info = -3;
764	return;
765	}
766	a[offs+j,offs+j] = 1/a[offs+j,offs+j];
767	ajj = -a[offs+j,offs+j];
768	}
769	else
770	{
771	ajj = -1;
772	}
773
774	//
775	// Compute elements 1:j-1 of j-th column.
776	//
777	if( j>0 )
778	{
779	i1_ = (offs+0) - (0);
780	for(i_=0; i_<=j-1;i_++)
781	{
782	tmp[i_] = a[i_+i1_,offs+j];
783	}
784	for(i=0; i<=j-1; i++)
785	{
786	if( i<j-1 )
787	{
788	i1_ = (i+1)-(offs+i+1);
789	v = 0.0;
790	for(i_=offs+i+1; i_<=offs+j-1;i_++)
791	{
792	v += a[offs+i,i_]*tmp[i_+i1_];
793	}
794	}
795	else
796	{
797	v = 0;
798	}
799	if( !isunit )
800	{
801	a[offs+i,offs+j] = v+a[offs+i,offs+i]*tmp[i];
802	}
803	else
804	{
805	a[offs+i,offs+j] = v+tmp[i];
806	}
807	}
808	for(i_=offs+0; i_<=offs+j-1;i_++)
809	{
810	a[i_,offs+j] = ajj*a[i_,offs+j];
811	}
812	}
813	}
814	}
815	else
816	{
817
818	//
819	// Compute inverse of lower triangular matrix.
820	//
821	for(j=n-1; j>=0; j--)
822	{
823	if( !isunit )
824	{
825	if( (double)(a[offs+j,offs+j])==(double)(0) )
826	{
827	info = -3;
828	return;
829	}
830	a[offs+j,offs+j] = 1/a[offs+j,offs+j];
831	ajj = -a[offs+j,offs+j];
832	}
833	else
834	{
835	ajj = -1;
836	}
837	if( j<n-1 )
838	{
839
840	//
841	// Compute elements j+1:n of j-th column.
842	//
843	i1_ = (offs+j+1) - (j+1);
844	for(i_=j+1; i_<=n-1;i_++)
845	{
846	tmp[i_] = a[i_+i1_,offs+j];
847	}
848	for(i=j+1; i<=n-1; i++)
849	{
850	if( i>j+1 )
851	{
852	i1_ = (j+1)-(offs+j+1);
853	v = 0.0;
854	for(i_=offs+j+1; i_<=offs+i-1;i_++)
855	{
856	v += a[offs+i,i_]*tmp[i_+i1_];
857	}
858	}
859	else
860	{
861	v = 0;
862	}
863	if( !isunit )
864	{
865	a[offs+i,offs+j] = v+a[offs+i,offs+i]*tmp[i];
866	}
867	else
868	{
869	a[offs+i,offs+j] = v+tmp[i];
870	}
871	}
872	for(i_=offs+j+1; i_<=offs+n-1;i_++)
873	{
874	a[i_,offs+j] = ajj*a[i_,offs+j];
875	}
876	}
877	}
878	}
879	return;
880	}
881
882	//
883	// Recursive case
884	//
885	ablas.ablassplitlength(ref a, n, ref n1, ref n2);
886	if( n2>0 )
887	{
888	if( isupper )
889	{
890	for(i=0; i<=n1-1; i++)
891	{
892	for(i_=offs+n1; i_<=offs+n-1;i_++)
893	{
894	a[offs+i,i_] = -1*a[offs+i,i_];
895	}
896	}
897	ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, isunit, 0, ref a, offs, offs+n1);
898	ablas.rmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, isunit, 0, ref a, offs, offs+n1);
899	}
900	else
901	{
902	for(i=0; i<=n2-1; i++)
903	{
904	for(i_=offs; i_<=offs+n1-1;i_++)
905	{
906	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
907	}
908	}
909	ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, isunit, 0, ref a, offs+n1, offs);
910	ablas.rmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, isunit, 0, ref a, offs+n1, offs);
911	}
912	rmatrixtrinverserec(ref a, offs+n1, n2, isupper, isunit, ref tmp, ref info, ref rep);
913	}
914	rmatrixtrinverserec(ref a, offs, n1, isupper, isunit, ref tmp, ref info, ref rep);
915	}
916
917
918	/*************************************************************************
919	Triangular matrix inversion, recursive subroutine
920
921	-- ALGLIB --
922	05.02.2010, Bochkanov Sergey.
923	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
924	Courant Institute, Argonne National Lab, and Rice University
925	February 29, 1992.
926	*************************************************************************/
927	private static void cmatrixtrinverserec(ref AP.Complex[,] a,
928	int offs,
929	int n,
930	bool isupper,
931	bool isunit,
932	ref AP.Complex[] tmp,
933	ref int info,
934	ref matinvreport rep)
935	{
936	int n1 = 0;
937	int n2 = 0;
938	int i = 0;
939	int j = 0;
940	AP.Complex v = 0;
941	AP.Complex ajj = 0;
942	int i_ = 0;
943	int i1_ = 0;
944
945	if( n<1 )
946	{
947	info = -1;
948	return;
949	}
950
951	//
952	// Base case
953	//
954	if( n<=ablas.ablascomplexblocksize(ref a) )
955	{
956	if( isupper )
957	{
958
959	//
960	// Compute inverse of upper triangular matrix.
961	//
962	for(j=0; j<=n-1; j++)
963	{
964	if( !isunit )
965	{
966	if( a[offs+j,offs+j]==0 )
967	{
968	info = -3;
969	return;
970	}
971	a[offs+j,offs+j] = 1/a[offs+j,offs+j];
972	ajj = -a[offs+j,offs+j];
973	}
974	else
975	{
976	ajj = -1;
977	}
978
979	//
980	// Compute elements 1:j-1 of j-th column.
981	//
982	if( j>0 )
983	{
984	i1_ = (offs+0) - (0);
985	for(i_=0; i_<=j-1;i_++)
986	{
987	tmp[i_] = a[i_+i1_,offs+j];
988	}
989	for(i=0; i<=j-1; i++)
990	{
991	if( i<j-1 )
992	{
993	i1_ = (i+1)-(offs+i+1);
994	v = 0.0;
995	for(i_=offs+i+1; i_<=offs+j-1;i_++)
996	{
997	v += a[offs+i,i_]*tmp[i_+i1_];
998	}
999	}
1000	else
1001	{
1002	v = 0;
1003	}
1004	if( !isunit )
1005	{
1006	a[offs+i,offs+j] = v+a[offs+i,offs+i]*tmp[i];
1007	}
1008	else
1009	{
1010	a[offs+i,offs+j] = v+tmp[i];
1011	}
1012	}
1013	for(i_=offs+0; i_<=offs+j-1;i_++)
1014	{
1015	a[i_,offs+j] = ajj*a[i_,offs+j];
1016	}
1017	}
1018	}
1019	}
1020	else
1021	{
1022
1023	//
1024	// Compute inverse of lower triangular matrix.
1025	//
1026	for(j=n-1; j>=0; j--)
1027	{
1028	if( !isunit )
1029	{
1030	if( a[offs+j,offs+j]==0 )
1031	{
1032	info = -3;
1033	return;
1034	}
1035	a[offs+j,offs+j] = 1/a[offs+j,offs+j];
1036	ajj = -a[offs+j,offs+j];
1037	}
1038	else
1039	{
1040	ajj = -1;
1041	}
1042	if( j<n-1 )
1043	{
1044
1045	//
1046	// Compute elements j+1:n of j-th column.
1047	//
1048	i1_ = (offs+j+1) - (j+1);
1049	for(i_=j+1; i_<=n-1;i_++)
1050	{
1051	tmp[i_] = a[i_+i1_,offs+j];
1052	}
1053	for(i=j+1; i<=n-1; i++)
1054	{
1055	if( i>j+1 )
1056	{
1057	i1_ = (j+1)-(offs+j+1);
1058	v = 0.0;
1059	for(i_=offs+j+1; i_<=offs+i-1;i_++)
1060	{
1061	v += a[offs+i,i_]*tmp[i_+i1_];
1062	}
1063	}
1064	else
1065	{
1066	v = 0;
1067	}
1068	if( !isunit )
1069	{
1070	a[offs+i,offs+j] = v+a[offs+i,offs+i]*tmp[i];
1071	}
1072	else
1073	{
1074	a[offs+i,offs+j] = v+tmp[i];
1075	}
1076	}
1077	for(i_=offs+j+1; i_<=offs+n-1;i_++)
1078	{
1079	a[i_,offs+j] = ajj*a[i_,offs+j];
1080	}
1081	}
1082	}
1083	}
1084	return;
1085	}
1086
1087	//
1088	// Recursive case
1089	//
1090	ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1091	if( n2>0 )
1092	{
1093	if( isupper )
1094	{
1095	for(i=0; i<=n1-1; i++)
1096	{
1097	for(i_=offs+n1; i_<=offs+n-1;i_++)
1098	{
1099	a[offs+i,i_] = -1*a[offs+i,i_];
1100	}
1101	}
1102	ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, isunit, 0, ref a, offs, offs+n1);
1103	ablas.cmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, isunit, 0, ref a, offs, offs+n1);
1104	}
1105	else
1106	{
1107	for(i=0; i<=n2-1; i++)
1108	{
1109	for(i_=offs; i_<=offs+n1-1;i_++)
1110	{
1111	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
1112	}
1113	}
1114	ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, isunit, 0, ref a, offs+n1, offs);
1115	ablas.cmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, isunit, 0, ref a, offs+n1, offs);
1116	}
1117	cmatrixtrinverserec(ref a, offs+n1, n2, isupper, isunit, ref tmp, ref info, ref rep);
1118	}
1119	cmatrixtrinverserec(ref a, offs, n1, isupper, isunit, ref tmp, ref info, ref rep);
1120	}
1121
1122
1123	private static void rmatrixluinverserec(ref double[,] a,
1124	int offs,
1125	int n,
1126	ref double[] work,
1127	ref int info,
1128	ref matinvreport rep)
1129	{
1130	int i = 0;
1131	int iws = 0;
1132	int j = 0;
1133	int jb = 0;
1134	int jj = 0;
1135	int jp = 0;
1136	int k = 0;
1137	double v = 0;
1138	int n1 = 0;
1139	int n2 = 0;
1140	int i_ = 0;
1141	int i1_ = 0;
1142
1143	if( n<1 )
1144	{
1145	info = -1;
1146	return;
1147	}
1148
1149	//
1150	// Base case
1151	//
1152	if( n<=ablas.ablasblocksize(ref a) )
1153	{
1154
1155	//
1156	// Form inv(U)
1157	//
1158	rmatrixtrinverserec(ref a, offs, n, true, false, ref work, ref info, ref rep);
1159	if( info<=0 )
1160	{
1161	return;
1162	}
1163
1164	//
1165	// Solve the equation inv(A)*L = inv(U) for inv(A).
1166	//
1167	for(j=n-1; j>=0; j--)
1168	{
1169
1170	//
1171	// Copy current column of L to WORK and replace with zeros.
1172	//
1173	for(i=j+1; i<=n-1; i++)
1174	{
1175	work[i] = a[offs+i,offs+j];
1176	a[offs+i,offs+j] = 0;
1177	}
1178
1179	//
1180	// Compute current column of inv(A).
1181	//
1182	if( j<n-1 )
1183	{
1184	for(i=0; i<=n-1; i++)
1185	{
1186	i1_ = (j+1)-(offs+j+1);
1187	v = 0.0;
1188	for(i_=offs+j+1; i_<=offs+n-1;i_++)
1189	{
1190	v += a[offs+i,i_]*work[i_+i1_];
1191	}
1192	a[offs+i,offs+j] = a[offs+i,offs+j]-v;
1193	}
1194	}
1195	}
1196	return;
1197	}
1198
1199	//
1200	// Recursive code:
1201	//
1202	// ( L1 ) ( U1 U12 )
1203	// A = ( ) * ( )
1204	// ( L12 L2 ) ( U2 )
1205	//
1206	// ( W X )
1207	// A^-1 = ( )
1208	// ( Y Z )
1209	//
1210	ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1211	System.Diagnostics.Debug.Assert(n2>0, "LUInverseRec: internal error!");
1212
1213	//
1214	// X := inv(U1)U12inv(U2)
1215	//
1216	ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, true, false, 0, ref a, offs, offs+n1);
1217	ablas.rmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, true, false, 0, ref a, offs, offs+n1);
1218
1219	//
1220	// Y := inv(L2)L12inv(L1)
1221	//
1222	ablas.rmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, false, true, 0, ref a, offs+n1, offs);
1223	ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, false, true, 0, ref a, offs+n1, offs);
1224
1225	//
1226	// W := inv(L1U1)+XY
1227	//
1228	rmatrixluinverserec(ref a, offs, n1, ref work, ref info, ref rep);
1229	if( info<=0 )
1230	{
1231	return;
1232	}
1233	ablas.rmatrixgemm(n1, n1, n2, 1.0, ref a, offs, offs+n1, 0, ref a, offs+n1, offs, 0, 1.0, ref a, offs, offs);
1234
1235	//
1236	// X := -X*inv(L2)
1237	// Y := -inv(U2)*Y
1238	//
1239	ablas.rmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, false, true, 0, ref a, offs, offs+n1);
1240	for(i=0; i<=n1-1; i++)
1241	{
1242	for(i_=offs+n1; i_<=offs+n-1;i_++)
1243	{
1244	a[offs+i,i_] = -1*a[offs+i,i_];
1245	}
1246	}
1247	ablas.rmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, true, false, 0, ref a, offs+n1, offs);
1248	for(i=0; i<=n2-1; i++)
1249	{
1250	for(i_=offs; i_<=offs+n1-1;i_++)
1251	{
1252	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
1253	}
1254	}
1255
1256	//
1257	// Z := inv(L2*U2)
1258	//
1259	rmatrixluinverserec(ref a, offs+n1, n2, ref work, ref info, ref rep);
1260	}
1261
1262
1263	private static void cmatrixluinverserec(ref AP.Complex[,] a,
1264	int offs,
1265	int n,
1266	ref AP.Complex[] work,
1267	ref int info,
1268	ref matinvreport rep)
1269	{
1270	int i = 0;
1271	int iws = 0;
1272	int j = 0;
1273	int jb = 0;
1274	int jj = 0;
1275	int jp = 0;
1276	int k = 0;
1277	AP.Complex v = 0;
1278	int n1 = 0;
1279	int n2 = 0;
1280	int i_ = 0;
1281	int i1_ = 0;
1282
1283	if( n<1 )
1284	{
1285	info = -1;
1286	return;
1287	}
1288
1289	//
1290	// Base case
1291	//
1292	if( n<=ablas.ablascomplexblocksize(ref a) )
1293	{
1294
1295	//
1296	// Form inv(U)
1297	//
1298	cmatrixtrinverserec(ref a, offs, n, true, false, ref work, ref info, ref rep);
1299	if( info<=0 )
1300	{
1301	return;
1302	}
1303
1304	//
1305	// Solve the equation inv(A)*L = inv(U) for inv(A).
1306	//
1307	for(j=n-1; j>=0; j--)
1308	{
1309
1310	//
1311	// Copy current column of L to WORK and replace with zeros.
1312	//
1313	for(i=j+1; i<=n-1; i++)
1314	{
1315	work[i] = a[offs+i,offs+j];
1316	a[offs+i,offs+j] = 0;
1317	}
1318
1319	//
1320	// Compute current column of inv(A).
1321	//
1322	if( j<n-1 )
1323	{
1324	for(i=0; i<=n-1; i++)
1325	{
1326	i1_ = (j+1)-(offs+j+1);
1327	v = 0.0;
1328	for(i_=offs+j+1; i_<=offs+n-1;i_++)
1329	{
1330	v += a[offs+i,i_]*work[i_+i1_];
1331	}
1332	a[offs+i,offs+j] = a[offs+i,offs+j]-v;
1333	}
1334	}
1335	}
1336	return;
1337	}
1338
1339	//
1340	// Recursive code:
1341	//
1342	// ( L1 ) ( U1 U12 )
1343	// A = ( ) * ( )
1344	// ( L12 L2 ) ( U2 )
1345	//
1346	// ( W X )
1347	// A^-1 = ( )
1348	// ( Y Z )
1349	//
1350	ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1351	System.Diagnostics.Debug.Assert(n2>0, "LUInverseRec: internal error!");
1352
1353	//
1354	// X := inv(U1)U12inv(U2)
1355	//
1356	ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, true, false, 0, ref a, offs, offs+n1);
1357	ablas.cmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, true, false, 0, ref a, offs, offs+n1);
1358
1359	//
1360	// Y := inv(L2)L12inv(L1)
1361	//
1362	ablas.cmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, false, true, 0, ref a, offs+n1, offs);
1363	ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, false, true, 0, ref a, offs+n1, offs);
1364
1365	//
1366	// W := inv(L1U1)+XY
1367	//
1368	cmatrixluinverserec(ref a, offs, n1, ref work, ref info, ref rep);
1369	if( info<=0 )
1370	{
1371	return;
1372	}
1373	ablas.cmatrixgemm(n1, n1, n2, 1.0, ref a, offs, offs+n1, 0, ref a, offs+n1, offs, 0, 1.0, ref a, offs, offs);
1374
1375	//
1376	// X := -X*inv(L2)
1377	// Y := -inv(U2)*Y
1378	//
1379	ablas.cmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, false, true, 0, ref a, offs, offs+n1);
1380	for(i=0; i<=n1-1; i++)
1381	{
1382	for(i_=offs+n1; i_<=offs+n-1;i_++)
1383	{
1384	a[offs+i,i_] = -1*a[offs+i,i_];
1385	}
1386	}
1387	ablas.cmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, true, false, 0, ref a, offs+n1, offs);
1388	for(i=0; i<=n2-1; i++)
1389	{
1390	for(i_=offs; i_<=offs+n1-1;i_++)
1391	{
1392	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
1393	}
1394	}
1395
1396	//
1397	// Z := inv(L2*U2)
1398	//
1399	cmatrixluinverserec(ref a, offs+n1, n2, ref work, ref info, ref rep);
1400	}
1401
1402
1403	/*************************************************************************
1404	Recursive subroutine for SPD inversion.
1405
1406	-- ALGLIB routine --
1407	10.02.2010
1408	Bochkanov Sergey
1409	*************************************************************************/
1410	private static void spdmatrixcholeskyinverserec(ref double[,] a,
1411	int offs,
1412	int n,
1413	bool isupper,
1414	ref double[] tmp)
1415	{
1416	int i = 0;
1417	int j = 0;
1418	double v = 0;
1419	int n1 = 0;
1420	int n2 = 0;
1421	int info2 = 0;
1422	matinvreport rep2 = new matinvreport();
1423	int i_ = 0;
1424	int i1_ = 0;
1425
1426	if( n<1 )
1427	{
1428	return;
1429	}
1430
1431	//
1432	// Base case
1433	//
1434	if( n<=ablas.ablasblocksize(ref a) )
1435	{
1436	rmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
1437	if( isupper )
1438	{
1439
1440	//
1441	// Compute the product U * U'.
1442	// NOTE: we never assume that diagonal of U is real
1443	//
1444	for(i=0; i<=n-1; i++)
1445	{
1446	if( i==0 )
1447	{
1448
1449	//
1450	// 1x1 matrix
1451	//
1452	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i]);
1453	}
1454	else
1455	{
1456
1457	//
1458	// (I+1)x(I+1) matrix,
1459	//
1460	// ( A11 A12 ) ( A11^H ) ( A11A11^H+A12A12^H A12*A22^H )
1461	// ( ) * ( ) = ( )
1462	// ( A22 ) ( A12^H A22^H ) ( A22A12^H A22A22^H )
1463	//
1464	// A11 is IxI, A22 is 1x1.
1465	//
1466	i1_ = (offs) - (0);
1467	for(i_=0; i_<=i-1;i_++)
1468	{
1469	tmp[i_] = a[i_+i1_,offs+i];
1470	}
1471	for(j=0; j<=i-1; j++)
1472	{
1473	v = a[offs+j,offs+i];
1474	i1_ = (j) - (offs+j);
1475	for(i_=offs+j; i_<=offs+i-1;i_++)
1476	{
1477	a[offs+j,i_] = a[offs+j,i_] + v*tmp[i_+i1_];
1478	}
1479	}
1480	v = a[offs+i,offs+i];
1481	for(i_=offs; i_<=offs+i-1;i_++)
1482	{
1483	a[i_,offs+i] = v*a[i_,offs+i];
1484	}
1485	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i]);
1486	}
1487	}
1488	}
1489	else
1490	{
1491
1492	//
1493	// Compute the product L' * L
1494	// NOTE: we never assume that diagonal of L is real
1495	//
1496	for(i=0; i<=n-1; i++)
1497	{
1498	if( i==0 )
1499	{
1500
1501	//
1502	// 1x1 matrix
1503	//
1504	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i]);
1505	}
1506	else
1507	{
1508
1509	//
1510	// (I+1)x(I+1) matrix,
1511	//
1512	// ( A11^H A21^H ) ( A11 ) ( A11^HA11+A21^HA21 A21^H*A22 )
1513	// ( ) * ( ) = ( )
1514	// ( A22^H ) ( A21 A22 ) ( A22^HA21 A22^HA22 )
1515	//
1516	// A11 is IxI, A22 is 1x1.
1517	//
1518	i1_ = (offs) - (0);
1519	for(i_=0; i_<=i-1;i_++)
1520	{
1521	tmp[i_] = a[offs+i,i_+i1_];
1522	}
1523	for(j=0; j<=i-1; j++)
1524	{
1525	v = a[offs+i,offs+j];
1526	i1_ = (0) - (offs);
1527	for(i_=offs; i_<=offs+j;i_++)
1528	{
1529	a[offs+j,i_] = a[offs+j,i_] + v*tmp[i_+i1_];
1530	}
1531	}
1532	v = a[offs+i,offs+i];
1533	for(i_=offs; i_<=offs+i-1;i_++)
1534	{
1535	a[offs+i,i_] = v*a[offs+i,i_];
1536	}
1537	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i]);
1538	}
1539	}
1540	}
1541	return;
1542	}
1543
1544	//
1545	// Recursive code: triangular factor inversion merged with
1546	// UU' or L'L multiplication
1547	//
1548	ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1549
1550	//
1551	// form off-diagonal block of trangular inverse
1552	//
1553	if( isupper )
1554	{
1555	for(i=0; i<=n1-1; i++)
1556	{
1557	for(i_=offs+n1; i_<=offs+n-1;i_++)
1558	{
1559	a[offs+i,i_] = -1*a[offs+i,i_];
1560	}
1561	}
1562	ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 0, ref a, offs, offs+n1);
1563	ablas.rmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, false, 0, ref a, offs, offs+n1);
1564	}
1565	else
1566	{
1567	for(i=0; i<=n2-1; i++)
1568	{
1569	for(i_=offs; i_<=offs+n1-1;i_++)
1570	{
1571	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
1572	}
1573	}
1574	ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 0, ref a, offs+n1, offs);
1575	ablas.rmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, false, 0, ref a, offs+n1, offs);
1576	}
1577
1578	//
1579	// invert first diagonal block
1580	//
1581	spdmatrixcholeskyinverserec(ref a, offs, n1, isupper, ref tmp);
1582
1583	//
1584	// update first diagonal block with off-diagonal block,
1585	// update off-diagonal block
1586	//
1587	if( isupper )
1588	{
1589	ablas.rmatrixsyrk(n1, n2, 1.0, ref a, offs, offs+n1, 0, 1.0, ref a, offs, offs, isupper);
1590	ablas.rmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, false, 1, ref a, offs, offs+n1);
1591	}
1592	else
1593	{
1594	ablas.rmatrixsyrk(n1, n2, 1.0, ref a, offs+n1, offs, 1, 1.0, ref a, offs, offs, isupper);
1595	ablas.rmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, false, 1, ref a, offs+n1, offs);
1596	}
1597
1598	//
1599	// invert second diagonal block
1600	//
1601	spdmatrixcholeskyinverserec(ref a, offs+n1, n2, isupper, ref tmp);
1602	}
1603
1604
1605	/*************************************************************************
1606	Recursive subroutine for HPD inversion.
1607
1608	-- ALGLIB routine --
1609	10.02.2010
1610	Bochkanov Sergey
1611	*************************************************************************/
1612	private static void hpdmatrixcholeskyinverserec(ref AP.Complex[,] a,
1613	int offs,
1614	int n,
1615	bool isupper,
1616	ref AP.Complex[] tmp)
1617	{
1618	int i = 0;
1619	int j = 0;
1620	AP.Complex v = 0;
1621	int n1 = 0;
1622	int n2 = 0;
1623	int info2 = 0;
1624	matinvreport rep2 = new matinvreport();
1625	int i_ = 0;
1626	int i1_ = 0;
1627
1628	if( n<1 )
1629	{
1630	return;
1631	}
1632
1633	//
1634	// Base case
1635	//
1636	if( n<=ablas.ablascomplexblocksize(ref a) )
1637	{
1638	cmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
1639	if( isupper )
1640	{
1641
1642	//
1643	// Compute the product U * U'.
1644	// NOTE: we never assume that diagonal of U is real
1645	//
1646	for(i=0; i<=n-1; i++)
1647	{
1648	if( i==0 )
1649	{
1650
1651	//
1652	// 1x1 matrix
1653	//
1654	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i].x)+AP.Math.Sqr(a[offs+i,offs+i].y);
1655	}
1656	else
1657	{
1658
1659	//
1660	// (I+1)x(I+1) matrix,
1661	//
1662	// ( A11 A12 ) ( A11^H ) ( A11A11^H+A12A12^H A12*A22^H )
1663	// ( ) * ( ) = ( )
1664	// ( A22 ) ( A12^H A22^H ) ( A22A12^H A22A22^H )
1665	//
1666	// A11 is IxI, A22 is 1x1.
1667	//
1668	i1_ = (offs) - (0);
1669	for(i_=0; i_<=i-1;i_++)
1670	{
1671	tmp[i_] = AP.Math.Conj(a[i_+i1_,offs+i]);
1672	}
1673	for(j=0; j<=i-1; j++)
1674	{
1675	v = a[offs+j,offs+i];
1676	i1_ = (j) - (offs+j);
1677	for(i_=offs+j; i_<=offs+i-1;i_++)
1678	{
1679	a[offs+j,i_] = a[offs+j,i_] + v*tmp[i_+i1_];
1680	}
1681	}
1682	v = AP.Math.Conj(a[offs+i,offs+i]);
1683	for(i_=offs; i_<=offs+i-1;i_++)
1684	{
1685	a[i_,offs+i] = v*a[i_,offs+i];
1686	}
1687	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i].x)+AP.Math.Sqr(a[offs+i,offs+i].y);
1688	}
1689	}
1690	}
1691	else
1692	{
1693
1694	//
1695	// Compute the product L' * L
1696	// NOTE: we never assume that diagonal of L is real
1697	//
1698	for(i=0; i<=n-1; i++)
1699	{
1700	if( i==0 )
1701	{
1702
1703	//
1704	// 1x1 matrix
1705	//
1706	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i].x)+AP.Math.Sqr(a[offs+i,offs+i].y);
1707	}
1708	else
1709	{
1710
1711	//
1712	// (I+1)x(I+1) matrix,
1713	//
1714	// ( A11^H A21^H ) ( A11 ) ( A11^HA11+A21^HA21 A21^H*A22 )
1715	// ( ) * ( ) = ( )
1716	// ( A22^H ) ( A21 A22 ) ( A22^HA21 A22^HA22 )
1717	//
1718	// A11 is IxI, A22 is 1x1.
1719	//
1720	i1_ = (offs) - (0);
1721	for(i_=0; i_<=i-1;i_++)
1722	{
1723	tmp[i_] = a[offs+i,i_+i1_];
1724	}
1725	for(j=0; j<=i-1; j++)
1726	{
1727	v = AP.Math.Conj(a[offs+i,offs+j]);
1728	i1_ = (0) - (offs);
1729	for(i_=offs; i_<=offs+j;i_++)
1730	{
1731	a[offs+j,i_] = a[offs+j,i_] + v*tmp[i_+i1_];
1732	}
1733	}
1734	v = AP.Math.Conj(a[offs+i,offs+i]);
1735	for(i_=offs; i_<=offs+i-1;i_++)
1736	{
1737	a[offs+i,i_] = v*a[offs+i,i_];
1738	}
1739	a[offs+i,offs+i] = AP.Math.Sqr(a[offs+i,offs+i].x)+AP.Math.Sqr(a[offs+i,offs+i].y);
1740	}
1741	}
1742	}
1743	return;
1744	}
1745
1746	//
1747	// Recursive code: triangular factor inversion merged with
1748	// UU' or L'L multiplication
1749	//
1750	ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1751
1752	//
1753	// form off-diagonal block of trangular inverse
1754	//
1755	if( isupper )
1756	{
1757	for(i=0; i<=n1-1; i++)
1758	{
1759	for(i_=offs+n1; i_<=offs+n-1;i_++)
1760	{
1761	a[offs+i,i_] = -1*a[offs+i,i_];
1762	}
1763	}
1764	ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 0, ref a, offs, offs+n1);
1765	ablas.cmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, false, 0, ref a, offs, offs+n1);
1766	}
1767	else
1768	{
1769	for(i=0; i<=n2-1; i++)
1770	{
1771	for(i_=offs; i_<=offs+n1-1;i_++)
1772	{
1773	a[offs+n1+i,i_] = -1*a[offs+n1+i,i_];
1774	}
1775	}
1776	ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 0, ref a, offs+n1, offs);
1777	ablas.cmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, false, 0, ref a, offs+n1, offs);
1778	}
1779
1780	//
1781	// invert first diagonal block
1782	//
1783	hpdmatrixcholeskyinverserec(ref a, offs, n1, isupper, ref tmp);
1784
1785	//
1786	// update first diagonal block with off-diagonal block,
1787	// update off-diagonal block
1788	//
1789	if( isupper )
1790	{
1791	ablas.cmatrixsyrk(n1, n2, 1.0, ref a, offs, offs+n1, 0, 1.0, ref a, offs, offs, isupper);
1792	ablas.cmatrixrighttrsm(n1, n2, ref a, offs+n1, offs+n1, isupper, false, 2, ref a, offs, offs+n1);
1793	}
1794	else
1795	{
1796	ablas.cmatrixsyrk(n1, n2, 1.0, ref a, offs+n1, offs, 2, 1.0, ref a, offs, offs, isupper);
1797	ablas.cmatrixlefttrsm(n2, n1, ref a, offs+n1, offs+n1, isupper, false, 2, ref a, offs+n1, offs);
1798	}
1799
1800	//
1801	// invert second diagonal block
1802	//
1803	hpdmatrixcholeskyinverserec(ref a, offs+n1, n2, isupper, ref tmp);
1804	}
1805	}
1806	}

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

Download in other formats:

Update cookies preferences