Context Navigation

source: trunk/sources/ALGLIB/sinverse.cs @ 2567

Visit:

Last change on this file since 2567 was 2563, checked in by gkronber, 15 years ago
Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 34.5 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 sinverse
32	{
33	/*************************************************************************
34	Inversion of a symmetric indefinite matrix
35
36	The algorithm gets an LDLT-decomposition as an input, generates matrix A^-1
37	and saves the lower or upper triangle of an inverse matrix depending on the
38	input (UDU' or LDL').
39
40	Input parameters:
41	A - LDLT-decomposition of the matrix,
42	Output of subroutine SMatrixLDLT.
43	N - size of matrix A.
44	IsUpper - storage format. If IsUpper = True, then the symmetric matrix
45	is given as decomposition A = UDU' and this decomposition
46	is stored in the upper triangle of matrix A and on the main
47	diagonal, and the lower triangle of matrix A is not used.
48	Pivots - a table of permutations, output of subroutine SMatrixLDLT.
49
50	Output parameters:
51	A - inverse of the matrix, whose LDLT-decomposition was stored
52	in matrix A as a subroutine input.
53	Array with elements [0..N-1, 0..N-1].
54	If IsUpper = True, then A contains the upper triangle of
55	matrix A^-1, and the elements below the main diagonal are
56	not used nor changed. The same applies if IsUpper = False.
57
58	Result:
59	True, if the matrix is not singular.
60	False, if the matrix is singular and could not be inverted.
61
62	-- LAPACK routine (version 3.0) --
63	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
64	Courant Institute, Argonne National Lab, and Rice University
65	March 31, 1993
66	*************************************************************************/
67	public static bool smatrixldltinverse(ref double[,] a,
68	ref int[] pivots,
69	int n,
70	bool isupper)
71	{
72	bool result = new bool();
73	double[] work = new double[0];
74	double[] work2 = new double[0];
75	int i = 0;
76	int k = 0;
77	int kp = 0;
78	int kstep = 0;
79	double ak = 0;
80	double akkp1 = 0;
81	double akp1 = 0;
82	double d = 0;
83	double t = 0;
84	double temp = 0;
85	int km1 = 0;
86	int kp1 = 0;
87	int l = 0;
88	int i1 = 0;
89	int i2 = 0;
90	double v = 0;
91	int i_ = 0;
92	int i1_ = 0;
93
94	work = new double[n+1];
95	work2 = new double[n+1];
96	result = true;
97
98	//
99	// Quick return if possible
100	//
101	if( n==0 )
102	{
103	return result;
104	}
105
106	//
107	// Check that the diagonal matrix D is nonsingular.
108	//
109	for(i=0; i<=n-1; i++)
110	{
111	if( pivots[i]>=0 & (double)(a[i,i])==(double)(0) )
112	{
113	result = false;
114	return result;
115	}
116	}
117	if( isupper )
118	{
119
120	//
121	// Compute inv(A) from the factorization A = UDU'.
122	//
123	// K+1 is the main loop index, increasing from 1 to N in steps of
124	// 1 or 2, depending on the size of the diagonal blocks.
125	//
126	k = 0;
127	while( k<=n-1 )
128	{
129	if( pivots[k]>=0 )
130	{
131
132	//
133	// 1 x 1 diagonal block
134	//
135	// Invert the diagonal block.
136	//
137	a[k,k] = 1/a[k,k];
138
139	//
140	// Compute column K+1 of the inverse.
141	//
142	if( k>0 )
143	{
144	i1_ = (0) - (1);
145	for(i_=1; i_<=k;i_++)
146	{
147	work[i_] = a[i_+i1_,k];
148	}
149	sblas.symmetricmatrixvectormultiply(ref a, isupper, 1-1, k+1-1-1, ref work, -1, ref work2);
150	i1_ = (1) - (0);
151	for(i_=0; i_<=k-1;i_++)
152	{
153	a[i_,k] = work2[i_+i1_];
154	}
155	v = 0.0;
156	for(i_=1; i_<=k;i_++)
157	{
158	v += work2[i_]*work[i_];
159	}
160	a[k,k] = a[k,k]-v;
161	}
162	kstep = 1;
163	}
164	else
165	{
166
167	//
168	// 2 x 2 diagonal block
169	//
170	// Invert the diagonal block.
171	//
172	t = Math.Abs(a[k,k+1]);
173	ak = a[k,k]/t;
174	akp1 = a[k+1,k+1]/t;
175	akkp1 = a[k,k+1]/t;
176	d = t(akakp1-1);
177	a[k,k] = akp1/d;
178	a[k+1,k+1] = ak/d;
179	a[k,k+1] = -(akkp1/d);
180
181	//
182	// Compute columns K+1 and K+1+1 of the inverse.
183	//
184	if( k>0 )
185	{
186	i1_ = (0) - (1);
187	for(i_=1; i_<=k;i_++)
188	{
189	work[i_] = a[i_+i1_,k];
190	}
191	sblas.symmetricmatrixvectormultiply(ref a, isupper, 0, k-1, ref work, -1, ref work2);
192	i1_ = (1) - (0);
193	for(i_=0; i_<=k-1;i_++)
194	{
195	a[i_,k] = work2[i_+i1_];
196	}
197	v = 0.0;
198	for(i_=1; i_<=k;i_++)
199	{
200	v += work[i_]*work2[i_];
201	}
202	a[k,k] = a[k,k]-v;
203	v = 0.0;
204	for(i_=0; i_<=k-1;i_++)
205	{
206	v += a[i_,k]*a[i_,k+1];
207	}
208	a[k,k+1] = a[k,k+1]-v;
209	i1_ = (0) - (1);
210	for(i_=1; i_<=k;i_++)
211	{
212	work[i_] = a[i_+i1_,k+1];
213	}
214	sblas.symmetricmatrixvectormultiply(ref a, isupper, 0, k-1, ref work, -1, ref work2);
215	i1_ = (1) - (0);
216	for(i_=0; i_<=k-1;i_++)
217	{
218	a[i_,k+1] = work2[i_+i1_];
219	}
220	v = 0.0;
221	for(i_=1; i_<=k;i_++)
222	{
223	v += work[i_]*work2[i_];
224	}
225	a[k+1,k+1] = a[k+1,k+1]-v;
226	}
227	kstep = 2;
228	}
229	if( pivots[k]>=0 )
230	{
231	kp = pivots[k];
232	}
233	else
234	{
235	kp = n+pivots[k];
236	}
237	if( kp!=k )
238	{
239
240	//
241	// Interchange rows and columns K and KP in the leading
242	// submatrix
243	//
244	i1_ = (0) - (1);
245	for(i_=1; i_<=kp;i_++)
246	{
247	work[i_] = a[i_+i1_,k];
248	}
249	for(i_=0; i_<=kp-1;i_++)
250	{
251	a[i_,k] = a[i_,kp];
252	}
253	i1_ = (1) - (0);
254	for(i_=0; i_<=kp-1;i_++)
255	{
256	a[i_,kp] = work[i_+i1_];
257	}
258	i1_ = (kp+1) - (1);
259	for(i_=1; i_<=k-1-kp;i_++)
260	{
261	work[i_] = a[i_+i1_,k];
262	}
263	for(i_=kp+1; i_<=k-1;i_++)
264	{
265	a[i_,k] = a[kp,i_];
266	}
267	i1_ = (1) - (kp+1);
268	for(i_=kp+1; i_<=k-1;i_++)
269	{
270	a[kp,i_] = work[i_+i1_];
271	}
272	temp = a[k,k];
273	a[k,k] = a[kp,kp];
274	a[kp,kp] = temp;
275	if( kstep==2 )
276	{
277	temp = a[k,k+1];
278	a[k,k+1] = a[kp,k+1];
279	a[kp,k+1] = temp;
280	}
281	}
282	k = k+kstep;
283	}
284	}
285	else
286	{
287
288	//
289	// Compute inv(A) from the factorization A = LDL'.
290	//
291	// K is the main loop index, increasing from 0 to N-1 in steps of
292	// 1 or 2, depending on the size of the diagonal blocks.
293	//
294	k = n-1;
295	while( k>=0 )
296	{
297	if( pivots[k]>=0 )
298	{
299
300	//
301	// 1 x 1 diagonal block
302	//
303	// Invert the diagonal block.
304	//
305	a[k,k] = 1/a[k,k];
306
307	//
308	// Compute column K+1 of the inverse.
309	//
310	if( k<n-1 )
311	{
312	i1_ = (k+1) - (1);
313	for(i_=1; i_<=n-k-1;i_++)
314	{
315	work[i_] = a[i_+i1_,k];
316	}
317	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n-1, ref work, -1, ref work2);
318	i1_ = (1) - (k+1);
319	for(i_=k+1; i_<=n-1;i_++)
320	{
321	a[i_,k] = work2[i_+i1_];
322	}
323	v = 0.0;
324	for(i_=1; i_<=n-k-1;i_++)
325	{
326	v += work[i_]*work2[i_];
327	}
328	a[k,k] = a[k,k]-v;
329	}
330	kstep = 1;
331	}
332	else
333	{
334
335	//
336	// 2 x 2 diagonal block
337	//
338	// Invert the diagonal block.
339	//
340	t = Math.Abs(a[k,k-1]);
341	ak = a[k-1,k-1]/t;
342	akp1 = a[k,k]/t;
343	akkp1 = a[k,k-1]/t;
344	d = t(akakp1-1);
345	a[k-1,k-1] = akp1/d;
346	a[k,k] = ak/d;
347	a[k,k-1] = -(akkp1/d);
348
349	//
350	// Compute columns K+1-1 and K+1 of the inverse.
351	//
352	if( k<n-1 )
353	{
354	i1_ = (k+1) - (1);
355	for(i_=1; i_<=n-k-1;i_++)
356	{
357	work[i_] = a[i_+i1_,k];
358	}
359	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n-1, ref work, -1, ref work2);
360	i1_ = (1) - (k+1);
361	for(i_=k+1; i_<=n-1;i_++)
362	{
363	a[i_,k] = work2[i_+i1_];
364	}
365	v = 0.0;
366	for(i_=1; i_<=n-k-1;i_++)
367	{
368	v += work[i_]*work2[i_];
369	}
370	a[k,k] = a[k,k]-v;
371	v = 0.0;
372	for(i_=k+1; i_<=n-1;i_++)
373	{
374	v += a[i_,k]*a[i_,k-1];
375	}
376	a[k,k-1] = a[k,k-1]-v;
377	i1_ = (k+1) - (1);
378	for(i_=1; i_<=n-k-1;i_++)
379	{
380	work[i_] = a[i_+i1_,k-1];
381	}
382	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n-1, ref work, -1, ref work2);
383	i1_ = (1) - (k+1);
384	for(i_=k+1; i_<=n-1;i_++)
385	{
386	a[i_,k-1] = work2[i_+i1_];
387	}
388	v = 0.0;
389	for(i_=1; i_<=n-k-1;i_++)
390	{
391	v += work[i_]*work2[i_];
392	}
393	a[k-1,k-1] = a[k-1,k-1]-v;
394	}
395	kstep = 2;
396	}
397	if( pivots[k]>=0 )
398	{
399	kp = pivots[k];
400	}
401	else
402	{
403	kp = pivots[k]+n;
404	}
405	if( kp!=k )
406	{
407
408	//
409	// Interchange rows and columns K and KP
410	//
411	if( kp<n-1 )
412	{
413	i1_ = (kp+1) - (1);
414	for(i_=1; i_<=n-kp-1;i_++)
415	{
416	work[i_] = a[i_+i1_,k];
417	}
418	for(i_=kp+1; i_<=n-1;i_++)
419	{
420	a[i_,k] = a[i_,kp];
421	}
422	i1_ = (1) - (kp+1);
423	for(i_=kp+1; i_<=n-1;i_++)
424	{
425	a[i_,kp] = work[i_+i1_];
426	}
427	}
428	i1_ = (k+1) - (1);
429	for(i_=1; i_<=kp-k-1;i_++)
430	{
431	work[i_] = a[i_+i1_,k];
432	}
433	for(i_=k+1; i_<=kp-1;i_++)
434	{
435	a[i_,k] = a[kp,i_];
436	}
437	i1_ = (1) - (k+1);
438	for(i_=k+1; i_<=kp-1;i_++)
439	{
440	a[kp,i_] = work[i_+i1_];
441	}
442	temp = a[k,k];
443	a[k,k] = a[kp,kp];
444	a[kp,kp] = temp;
445	if( kstep==2 )
446	{
447	temp = a[k,k-1];
448	a[k,k-1] = a[kp,k-1];
449	a[kp,k-1] = temp;
450	}
451	}
452	k = k-kstep;
453	}
454	}
455	return result;
456	}
457
458
459	/*************************************************************************
460	Inversion of a symmetric indefinite matrix
461
462	Given a lower or upper triangle of matrix A, the algorithm generates
463	matrix A^-1 and saves the lower or upper triangle depending on the input.
464
465	Input parameters:
466	A - matrix to be inverted (upper or lower triangle).
467	Array with elements [0..N-1, 0..N-1].
468	N - size of matrix A.
469	IsUpper - storage format. If IsUpper = True, then the upper
470	triangle of matrix A is given, otherwise the lower
471	triangle is given.
472
473	Output parameters:
474	A - inverse of matrix A.
475	Array with elements [0..N-1, 0..N-1].
476	If IsUpper = True, then A contains the upper triangle of
477	matrix A^-1, and the elements below the main diagonal are
478	not used nor changed.
479	The same applies if IsUpper = False.
480
481	Result:
482	True, if the matrix is not singular.
483	False, if the matrix is singular and could not be inverted.
484
485	-- LAPACK routine (version 3.0) --
486	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
487	Courant Institute, Argonne National Lab, and Rice University
488	March 31, 1993
489	*************************************************************************/
490	public static bool smatrixinverse(ref double[,] a,
491	int n,
492	bool isupper)
493	{
494	bool result = new bool();
495	int[] pivots = new int[0];
496
497	ldlt.smatrixldlt(ref a, n, isupper, ref pivots);
498	result = smatrixldltinverse(ref a, ref pivots, n, isupper);
499	return result;
500	}
501
502
503	public static bool inverseldlt(ref double[,] a,
504	ref int[] pivots,
505	int n,
506	bool isupper)
507	{
508	bool result = new bool();
509	double[] work = new double[0];
510	double[] work2 = new double[0];
511	int i = 0;
512	int k = 0;
513	int kp = 0;
514	int kstep = 0;
515	double ak = 0;
516	double akkp1 = 0;
517	double akp1 = 0;
518	double d = 0;
519	double t = 0;
520	double temp = 0;
521	int km1 = 0;
522	int kp1 = 0;
523	int l = 0;
524	int i1 = 0;
525	int i2 = 0;
526	double v = 0;
527	int i_ = 0;
528	int i1_ = 0;
529
530	work = new double[n+1];
531	work2 = new double[n+1];
532	result = true;
533
534	//
535	// Quick return if possible
536	//
537	if( n==0 )
538	{
539	return result;
540	}
541
542	//
543	// Check that the diagonal matrix D is nonsingular.
544	//
545	for(i=1; i<=n; i++)
546	{
547	if( pivots[i]>0 & (double)(a[i,i])==(double)(0) )
548	{
549	result = false;
550	return result;
551	}
552	}
553	if( isupper )
554	{
555
556	//
557	// Compute inv(A) from the factorization A = UDU'.
558	//
559	// K is the main loop index, increasing from 1 to N in steps of
560	// 1 or 2, depending on the size of the diagonal blocks.
561	//
562	k = 1;
563	while( k<=n )
564	{
565	if( pivots[k]>0 )
566	{
567
568	//
569	// 1 x 1 diagonal block
570	//
571	// Invert the diagonal block.
572	//
573	a[k,k] = 1/a[k,k];
574
575	//
576	// Compute column K of the inverse.
577	//
578	if( k>1 )
579	{
580	km1 = k-1;
581	for(i_=1; i_<=km1;i_++)
582	{
583	work[i_] = a[i_,k];
584	}
585	sblas.symmetricmatrixvectormultiply(ref a, isupper, 1, k-1, ref work, -1, ref work2);
586	for(i_=1; i_<=km1;i_++)
587	{
588	a[i_,k] = work2[i_];
589	}
590	v = 0.0;
591	for(i_=1; i_<=km1;i_++)
592	{
593	v += work2[i_]*work[i_];
594	}
595	a[k,k] = a[k,k]-v;
596	}
597	kstep = 1;
598	}
599	else
600	{
601
602	//
603	// 2 x 2 diagonal block
604	//
605	// Invert the diagonal block.
606	//
607	t = Math.Abs(a[k,k+1]);
608	ak = a[k,k]/t;
609	akp1 = a[k+1,k+1]/t;
610	akkp1 = a[k,k+1]/t;
611	d = t(akakp1-1);
612	a[k,k] = akp1/d;
613	a[k+1,k+1] = ak/d;
614	a[k,k+1] = -(akkp1/d);
615
616	//
617	// Compute columns K and K+1 of the inverse.
618	//
619	if( k>1 )
620	{
621	km1 = k-1;
622	kp1 = k+1;
623	for(i_=1; i_<=km1;i_++)
624	{
625	work[i_] = a[i_,k];
626	}
627	sblas.symmetricmatrixvectormultiply(ref a, isupper, 1, k-1, ref work, -1, ref work2);
628	for(i_=1; i_<=km1;i_++)
629	{
630	a[i_,k] = work2[i_];
631	}
632	v = 0.0;
633	for(i_=1; i_<=km1;i_++)
634	{
635	v += work[i_]*work2[i_];
636	}
637	a[k,k] = a[k,k]-v;
638	v = 0.0;
639	for(i_=1; i_<=km1;i_++)
640	{
641	v += a[i_,k]*a[i_,kp1];
642	}
643	a[k,k+1] = a[k,k+1]-v;
644	for(i_=1; i_<=km1;i_++)
645	{
646	work[i_] = a[i_,kp1];
647	}
648	sblas.symmetricmatrixvectormultiply(ref a, isupper, 1, k-1, ref work, -1, ref work2);
649	for(i_=1; i_<=km1;i_++)
650	{
651	a[i_,kp1] = work2[i_];
652	}
653	v = 0.0;
654	for(i_=1; i_<=km1;i_++)
655	{
656	v += work[i_]*work2[i_];
657	}
658	a[k+1,k+1] = a[k+1,k+1]-v;
659	}
660	kstep = 2;
661	}
662	kp = Math.Abs(pivots[k]);
663	if( kp!=k )
664	{
665
666	//
667	// Interchange rows and columns K and KP in the leading
668	// submatrix A(1:k+1,1:k+1)
669	//
670	l = kp-1;
671	for(i_=1; i_<=l;i_++)
672	{
673	work[i_] = a[i_,k];
674	}
675	for(i_=1; i_<=l;i_++)
676	{
677	a[i_,k] = a[i_,kp];
678	}
679	for(i_=1; i_<=l;i_++)
680	{
681	a[i_,kp] = work[i_];
682	}
683	l = k-kp-1;
684	i1 = kp+1;
685	i2 = k-1;
686	i1_ = (i1) - (1);
687	for(i_=1; i_<=l;i_++)
688	{
689	work[i_] = a[i_+i1_,k];
690	}
691	for(i_=i1; i_<=i2;i_++)
692	{
693	a[i_,k] = a[kp,i_];
694	}
695	i1_ = (1) - (i1);
696	for(i_=i1; i_<=i2;i_++)
697	{
698	a[kp,i_] = work[i_+i1_];
699	}
700	temp = a[k,k];
701	a[k,k] = a[kp,kp];
702	a[kp,kp] = temp;
703	if( kstep==2 )
704	{
705	temp = a[k,k+1];
706	a[k,k+1] = a[kp,k+1];
707	a[kp,k+1] = temp;
708	}
709	}
710	k = k+kstep;
711	}
712	}
713	else
714	{
715
716	//
717	// Compute inv(A) from the factorization A = LDL'.
718	//
719	// K is the main loop index, increasing from 1 to N in steps of
720	// 1 or 2, depending on the size of the diagonal blocks.
721	//
722	k = n;
723	while( k>=1 )
724	{
725	if( pivots[k]>0 )
726	{
727
728	//
729	// 1 x 1 diagonal block
730	//
731	// Invert the diagonal block.
732	//
733	a[k,k] = 1/a[k,k];
734
735	//
736	// Compute column K of the inverse.
737	//
738	if( k<n )
739	{
740	kp1 = k+1;
741	km1 = k-1;
742	l = n-k;
743	i1_ = (kp1) - (1);
744	for(i_=1; i_<=l;i_++)
745	{
746	work[i_] = a[i_+i1_,k];
747	}
748	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n, ref work, -1, ref work2);
749	i1_ = (1) - (kp1);
750	for(i_=kp1; i_<=n;i_++)
751	{
752	a[i_,k] = work2[i_+i1_];
753	}
754	v = 0.0;
755	for(i_=1; i_<=l;i_++)
756	{
757	v += work[i_]*work2[i_];
758	}
759	a[k,k] = a[k,k]-v;
760	}
761	kstep = 1;
762	}
763	else
764	{
765
766	//
767	// 2 x 2 diagonal block
768	//
769	// Invert the diagonal block.
770	//
771	t = Math.Abs(a[k,k-1]);
772	ak = a[k-1,k-1]/t;
773	akp1 = a[k,k]/t;
774	akkp1 = a[k,k-1]/t;
775	d = t(akakp1-1);
776	a[k-1,k-1] = akp1/d;
777	a[k,k] = ak/d;
778	a[k,k-1] = -(akkp1/d);
779
780	//
781	// Compute columns K-1 and K of the inverse.
782	//
783	if( k<n )
784	{
785	kp1 = k+1;
786	km1 = k-1;
787	l = n-k;
788	i1_ = (kp1) - (1);
789	for(i_=1; i_<=l;i_++)
790	{
791	work[i_] = a[i_+i1_,k];
792	}
793	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n, ref work, -1, ref work2);
794	i1_ = (1) - (kp1);
795	for(i_=kp1; i_<=n;i_++)
796	{
797	a[i_,k] = work2[i_+i1_];
798	}
799	v = 0.0;
800	for(i_=1; i_<=l;i_++)
801	{
802	v += work[i_]*work2[i_];
803	}
804	a[k,k] = a[k,k]-v;
805	v = 0.0;
806	for(i_=kp1; i_<=n;i_++)
807	{
808	v += a[i_,k]*a[i_,km1];
809	}
810	a[k,k-1] = a[k,k-1]-v;
811	i1_ = (kp1) - (1);
812	for(i_=1; i_<=l;i_++)
813	{
814	work[i_] = a[i_+i1_,km1];
815	}
816	sblas.symmetricmatrixvectormultiply(ref a, isupper, k+1, n, ref work, -1, ref work2);
817	i1_ = (1) - (kp1);
818	for(i_=kp1; i_<=n;i_++)
819	{
820	a[i_,km1] = work2[i_+i1_];
821	}
822	v = 0.0;
823	for(i_=1; i_<=l;i_++)
824	{
825	v += work[i_]*work2[i_];
826	}
827	a[k-1,k-1] = a[k-1,k-1]-v;
828	}
829	kstep = 2;
830	}
831	kp = Math.Abs(pivots[k]);
832	if( kp!=k )
833	{
834
835	//
836	// Interchange rows and columns K and KP in the trailing
837	// submatrix A(k-1:n,k-1:n)
838	//
839	if( kp<n )
840	{
841	l = n-kp;
842	kp1 = kp+1;
843	i1_ = (kp1) - (1);
844	for(i_=1; i_<=l;i_++)
845	{
846	work[i_] = a[i_+i1_,k];
847	}
848	for(i_=kp1; i_<=n;i_++)
849	{
850	a[i_,k] = a[i_,kp];
851	}
852	i1_ = (1) - (kp1);
853	for(i_=kp1; i_<=n;i_++)
854	{
855	a[i_,kp] = work[i_+i1_];
856	}
857	}
858	l = kp-k-1;
859	i1 = k+1;
860	i2 = kp-1;
861	i1_ = (i1) - (1);
862	for(i_=1; i_<=l;i_++)
863	{
864	work[i_] = a[i_+i1_,k];
865	}
866	for(i_=i1; i_<=i2;i_++)
867	{
868	a[i_,k] = a[kp,i_];
869	}
870	i1_ = (1) - (i1);
871	for(i_=i1; i_<=i2;i_++)
872	{
873	a[kp,i_] = work[i_+i1_];
874	}
875	temp = a[k,k];
876	a[k,k] = a[kp,kp];
877	a[kp,kp] = temp;
878	if( kstep==2 )
879	{
880	temp = a[k,k-1];
881	a[k,k-1] = a[kp,k-1];
882	a[kp,k-1] = temp;
883	}
884	}
885	k = k-kstep;
886	}
887	}
888	return result;
889	}
890
891
892	public static bool inversesymmetricindefinite(ref double[,] a,
893	int n,
894	bool isupper)
895	{
896	bool result = new bool();
897	int[] pivots = new int[0];
898
899	ldlt.ldltdecomposition(ref a, n, isupper, ref pivots);
900	result = inverseldlt(ref a, ref pivots, n, isupper);
901	return result;
902	}
903	}
904	}

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

Download in other formats:

Update cookies preferences