Context Navigation

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/ablas.cs @ 17712

Visit:

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

Line
1	/*************************************************************************
2	Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
3
4	>>> SOURCE LICENSE >>>
5	This program is free software; you can redistribute it and/or modify
6	it under the terms of the GNU General Public License as published by
7	the Free Software Foundation (www.fsf.org); either version 2 of the
8	License, or (at your option) any later version.
9
10	This program is distributed in the hope that it will be useful,
11	but WITHOUT ANY WARRANTY; without even the implied warranty of
12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	GNU General Public License for more details.
14
15	A copy of the GNU General Public License is available at
16	http://www.fsf.org/licensing/licenses
17
18	>>> END OF LICENSE >>>
19	*************************************************************************/
20
21	using System;
22
23	namespace alglib
24	{
25	public class ablas
26	{
27	/*************************************************************************
28	Splits matrix length in two parts, left part should match ABLAS block size
29
30	INPUT PARAMETERS
31	A - real matrix, is passed to ensure that we didn't split
32	complex matrix using real splitting subroutine.
33	matrix itself is not changed.
34	N - length, N>0
35
36	OUTPUT PARAMETERS
37	N1 - length
38	N2 - length
39
40	N1+N2=N, N1>=N2, N2 may be zero
41
42	-- ALGLIB routine --
43	15.12.2009
44	Bochkanov Sergey
45	*************************************************************************/
46	public static void ablassplitlength(ref double[,] a,
47	int n,
48	ref int n1,
49	ref int n2)
50	{
51	if( n>ablasblocksize(ref a) )
52	{
53	ablasinternalsplitlength(n, ablasblocksize(ref a), ref n1, ref n2);
54	}
55	else
56	{
57	ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
58	}
59	}
60
61
62	/*************************************************************************
63	Complex ABLASSplitLength
64
65	-- ALGLIB routine --
66	15.12.2009
67	Bochkanov Sergey
68	*************************************************************************/
69	public static void ablascomplexsplitlength(ref AP.Complex[,] a,
70	int n,
71	ref int n1,
72	ref int n2)
73	{
74	if( n>ablascomplexblocksize(ref a) )
75	{
76	ablasinternalsplitlength(n, ablascomplexblocksize(ref a), ref n1, ref n2);
77	}
78	else
79	{
80	ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
81	}
82	}
83
84
85	/*************************************************************************
86	Returns block size - subdivision size where cache-oblivious soubroutines
87	switch to the optimized kernel.
88
89	INPUT PARAMETERS
90	A - real matrix, is passed to ensure that we didn't split
91	complex matrix using real splitting subroutine.
92	matrix itself is not changed.
93
94	-- ALGLIB routine --
95	15.12.2009
96	Bochkanov Sergey
97	*************************************************************************/
98	public static int ablasblocksize(ref double[,] a)
99	{
100	int result = 0;
101
102	result = 32;
103	return result;
104	}
105
106
107	/*************************************************************************
108	Block size for complex subroutines.
109
110	-- ALGLIB routine --
111	15.12.2009
112	Bochkanov Sergey
113	*************************************************************************/
114	public static int ablascomplexblocksize(ref AP.Complex[,] a)
115	{
116	int result = 0;
117
118	result = 24;
119	return result;
120	}
121
122
123	/*************************************************************************
124	Microblock size
125
126	-- ALGLIB routine --
127	15.12.2009
128	Bochkanov Sergey
129	*************************************************************************/
130	public static int ablasmicroblocksize()
131	{
132	int result = 0;
133
134	result = 8;
135	return result;
136	}
137
138
139	/*************************************************************************
140	Cache-oblivous complex "copy-and-transpose"
141
142	Input parameters:
143	M - number of rows
144	N - number of columns
145	A - source matrix, MxN submatrix is copied and transposed
146	IA - submatrix offset (row index)
147	JA - submatrix offset (column index)
148	A - destination matrix
149	IB - submatrix offset (row index)
150	JB - submatrix offset (column index)
151	*************************************************************************/
152	public static void cmatrixtranspose(int m,
153	int n,
154	ref AP.Complex[,] a,
155	int ia,
156	int ja,
157	ref AP.Complex[,] b,
158	int ib,
159	int jb)
160	{
161	int i = 0;
162	int s1 = 0;
163	int s2 = 0;
164	int i_ = 0;
165	int i1_ = 0;
166
167	if( m<=2ablascomplexblocksize(ref a) & n<=2ablascomplexblocksize(ref a) )
168	{
169
170	//
171	// base case
172	//
173	for(i=0; i<=m-1; i++)
174	{
175	i1_ = (ja) - (ib);
176	for(i_=ib; i_<=ib+n-1;i_++)
177	{
178	b[i_,jb+i] = a[ia+i,i_+i1_];
179	}
180	}
181	}
182	else
183	{
184
185	//
186	// Cache-oblivious recursion
187	//
188	if( m>n )
189	{
190	ablascomplexsplitlength(ref a, m, ref s1, ref s2);
191	cmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
192	cmatrixtranspose(s2, n, ref a, ia+s1, ja, ref b, ib, jb+s1);
193	}
194	else
195	{
196	ablascomplexsplitlength(ref a, n, ref s1, ref s2);
197	cmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
198	cmatrixtranspose(m, s2, ref a, ia, ja+s1, ref b, ib+s1, jb);
199	}
200	}
201	}
202
203
204	/*************************************************************************
205	Cache-oblivous real "copy-and-transpose"
206
207	Input parameters:
208	M - number of rows
209	N - number of columns
210	A - source matrix, MxN submatrix is copied and transposed
211	IA - submatrix offset (row index)
212	JA - submatrix offset (column index)
213	A - destination matrix
214	IB - submatrix offset (row index)
215	JB - submatrix offset (column index)
216	*************************************************************************/
217	public static void rmatrixtranspose(int m,
218	int n,
219	ref double[,] a,
220	int ia,
221	int ja,
222	ref double[,] b,
223	int ib,
224	int jb)
225	{
226	int i = 0;
227	int s1 = 0;
228	int s2 = 0;
229	int i_ = 0;
230	int i1_ = 0;
231
232	if( m<=2ablasblocksize(ref a) & n<=2ablasblocksize(ref a) )
233	{
234
235	//
236	// base case
237	//
238	for(i=0; i<=m-1; i++)
239	{
240	i1_ = (ja) - (ib);
241	for(i_=ib; i_<=ib+n-1;i_++)
242	{
243	b[i_,jb+i] = a[ia+i,i_+i1_];
244	}
245	}
246	}
247	else
248	{
249
250	//
251	// Cache-oblivious recursion
252	//
253	if( m>n )
254	{
255	ablassplitlength(ref a, m, ref s1, ref s2);
256	rmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
257	rmatrixtranspose(s2, n, ref a, ia+s1, ja, ref b, ib, jb+s1);
258	}
259	else
260	{
261	ablassplitlength(ref a, n, ref s1, ref s2);
262	rmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
263	rmatrixtranspose(m, s2, ref a, ia, ja+s1, ref b, ib+s1, jb);
264	}
265	}
266	}
267
268
269	/*************************************************************************
270	Rank-1 correction: A := A + u*v'
271
272	INPUT PARAMETERS:
273	M - number of rows
274	N - number of columns
275	A - target matrix, MxN submatrix is updated
276	IA - submatrix offset (row index)
277	JA - submatrix offset (column index)
278	U - vector #1
279	IU - subvector offset
280	V - vector #2
281	IV - subvector offset
282	*************************************************************************/
283	public static void cmatrixrank1(int m,
284	int n,
285	ref AP.Complex[,] a,
286	int ia,
287	int ja,
288	ref AP.Complex[] u,
289	int iu,
290	ref AP.Complex[] v,
291	int iv)
292	{
293	int i = 0;
294	AP.Complex s = 0;
295	int i_ = 0;
296	int i1_ = 0;
297
298	if( m==0 \| n==0 )
299	{
300	return;
301	}
302	if( ablasf.cmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv) )
303	{
304	return;
305	}
306	for(i=0; i<=m-1; i++)
307	{
308	s = u[iu+i];
309	i1_ = (iv) - (ja);
310	for(i_=ja; i_<=ja+n-1;i_++)
311	{
312	a[ia+i,i_] = a[ia+i,i_] + s*v[i_+i1_];
313	}
314	}
315	}
316
317
318	/*************************************************************************
319	Rank-1 correction: A := A + u*v'
320
321	INPUT PARAMETERS:
322	M - number of rows
323	N - number of columns
324	A - target matrix, MxN submatrix is updated
325	IA - submatrix offset (row index)
326	JA - submatrix offset (column index)
327	U - vector #1
328	IU - subvector offset
329	V - vector #2
330	IV - subvector offset
331	*************************************************************************/
332	public static void rmatrixrank1(int m,
333	int n,
334	ref double[,] a,
335	int ia,
336	int ja,
337	ref double[] u,
338	int iu,
339	ref double[] v,
340	int iv)
341	{
342	int i = 0;
343	double s = 0;
344	int i_ = 0;
345	int i1_ = 0;
346
347	if( m==0 \| n==0 )
348	{
349	return;
350	}
351	if( ablasf.rmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv) )
352	{
353	return;
354	}
355	for(i=0; i<=m-1; i++)
356	{
357	s = u[iu+i];
358	i1_ = (iv) - (ja);
359	for(i_=ja; i_<=ja+n-1;i_++)
360	{
361	a[ia+i,i_] = a[ia+i,i_] + s*v[i_+i1_];
362	}
363	}
364	}
365
366
367	/*************************************************************************
368	Matrix-vector product: y := op(A)*x
369
370	INPUT PARAMETERS:
371	M - number of rows of op(A)
372	M>=0
373	N - number of columns of op(A)
374	N>=0
375	A - target matrix
376	IA - submatrix offset (row index)
377	JA - submatrix offset (column index)
378	OpA - operation type:
379	* OpA=0 => op(A) = A
380	* OpA=1 => op(A) = A^T
381	* OpA=2 => op(A) = A^H
382	X - input vector
383	IX - subvector offset
384	IY - subvector offset
385
386	OUTPUT PARAMETERS:
387	Y - vector which stores result
388
389	if M=0, then subroutine does nothing.
390	if N=0, Y is filled by zeros.
391
392
393	-- ALGLIB routine --
394
395	28.01.2010
396	Bochkanov Sergey
397	*************************************************************************/
398	public static void cmatrixmv(int m,
399	int n,
400	ref AP.Complex[,] a,
401	int ia,
402	int ja,
403	int opa,
404	ref AP.Complex[] x,
405	int ix,
406	ref AP.Complex[] y,
407	int iy)
408	{
409	int i = 0;
410	AP.Complex v = 0;
411	int i_ = 0;
412	int i1_ = 0;
413
414	if( m==0 )
415	{
416	return;
417	}
418	if( n==0 )
419	{
420	for(i=0; i<=m-1; i++)
421	{
422	y[iy+i] = 0;
423	}
424	return;
425	}
426	if( ablasf.cmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy) )
427	{
428	return;
429	}
430	if( opa==0 )
431	{
432
433	//
434	// y = A*x
435	//
436	for(i=0; i<=m-1; i++)
437	{
438	i1_ = (ix)-(ja);
439	v = 0.0;
440	for(i_=ja; i_<=ja+n-1;i_++)
441	{
442	v += a[ia+i,i_]*x[i_+i1_];
443	}
444	y[iy+i] = v;
445	}
446	return;
447	}
448	if( opa==1 )
449	{
450
451	//
452	// y = A^T*x
453	//
454	for(i=0; i<=m-1; i++)
455	{
456	y[iy+i] = 0;
457	}
458	for(i=0; i<=n-1; i++)
459	{
460	v = x[ix+i];
461	i1_ = (ja) - (iy);
462	for(i_=iy; i_<=iy+m-1;i_++)
463	{
464	y[i_] = y[i_] + v*a[ia+i,i_+i1_];
465	}
466	}
467	return;
468	}
469	if( opa==2 )
470	{
471
472	//
473	// y = A^H*x
474	//
475	for(i=0; i<=m-1; i++)
476	{
477	y[iy+i] = 0;
478	}
479	for(i=0; i<=n-1; i++)
480	{
481	v = x[ix+i];
482	i1_ = (ja) - (iy);
483	for(i_=iy; i_<=iy+m-1;i_++)
484	{
485	y[i_] = y[i_] + v*AP.Math.Conj(a[ia+i,i_+i1_]);
486	}
487	}
488	return;
489	}
490	}
491
492
493	/*************************************************************************
494	Matrix-vector product: y := op(A)*x
495
496	INPUT PARAMETERS:
497	M - number of rows of op(A)
498	N - number of columns of op(A)
499	A - target matrix
500	IA - submatrix offset (row index)
501	JA - submatrix offset (column index)
502	OpA - operation type:
503	* OpA=0 => op(A) = A
504	* OpA=1 => op(A) = A^T
505	X - input vector
506	IX - subvector offset
507	IY - subvector offset
508
509	OUTPUT PARAMETERS:
510	Y - vector which stores result
511
512	if M=0, then subroutine does nothing.
513	if N=0, Y is filled by zeros.
514
515
516	-- ALGLIB routine --
517
518	28.01.2010
519	Bochkanov Sergey
520	*************************************************************************/
521	public static void rmatrixmv(int m,
522	int n,
523	ref double[,] a,
524	int ia,
525	int ja,
526	int opa,
527	ref double[] x,
528	int ix,
529	ref double[] y,
530	int iy)
531	{
532	int i = 0;
533	double v = 0;
534	int i_ = 0;
535	int i1_ = 0;
536
537	if( m==0 )
538	{
539	return;
540	}
541	if( n==0 )
542	{
543	for(i=0; i<=m-1; i++)
544	{
545	y[iy+i] = 0;
546	}
547	return;
548	}
549	if( ablasf.rmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy) )
550	{
551	return;
552	}
553	if( opa==0 )
554	{
555
556	//
557	// y = A*x
558	//
559	for(i=0; i<=m-1; i++)
560	{
561	i1_ = (ix)-(ja);
562	v = 0.0;
563	for(i_=ja; i_<=ja+n-1;i_++)
564	{
565	v += a[ia+i,i_]*x[i_+i1_];
566	}
567	y[iy+i] = v;
568	}
569	return;
570	}
571	if( opa==1 )
572	{
573
574	//
575	// y = A^T*x
576	//
577	for(i=0; i<=m-1; i++)
578	{
579	y[iy+i] = 0;
580	}
581	for(i=0; i<=n-1; i++)
582	{
583	v = x[ix+i];
584	i1_ = (ja) - (iy);
585	for(i_=iy; i_<=iy+m-1;i_++)
586	{
587	y[i_] = y[i_] + v*a[ia+i,i_+i1_];
588	}
589	}
590	return;
591	}
592	}
593
594
595	/*************************************************************************
596	This subroutine calculates X*op(A^-1) where:
597	* X is MxN general matrix
598	* A is NxN upper/lower triangular/unitriangular matrix
599	* "op" may be identity transformation, transposition, conjugate transposition
600
601	Multiplication result replaces X.
602	Cache-oblivious algorithm is used.
603
604	INPUT PARAMETERS
605	N - matrix size, N>=0
606	M - matrix size, N>=0
607	A - matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1]
608	I1 - submatrix offset
609	J1 - submatrix offset
610	IsUpper - whether matrix is upper triangular
611	IsUnit - whether matrix is unitriangular
612	OpType - transformation type:
613	* 0 - no transformation
614	* 1 - transposition
615	* 2 - conjugate transposition
616	C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
617	I2 - submatrix offset
618	J2 - submatrix offset
619
620	-- ALGLIB routine --
621	15.12.2009
622	Bochkanov Sergey
623	*************************************************************************/
624	public static void cmatrixrighttrsm(int m,
625	int n,
626	ref AP.Complex[,] a,
627	int i1,
628	int j1,
629	bool isupper,
630	bool isunit,
631	int optype,
632	ref AP.Complex[,] x,
633	int i2,
634	int j2)
635	{
636	int s1 = 0;
637	int s2 = 0;
638	int bs = 0;
639
640	bs = ablascomplexblocksize(ref a);
641	if( m<=bs & n<=bs )
642	{
643	cmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
644	return;
645	}
646	if( m>=n )
647	{
648
649	//
650	// Split X: XA = (X1 X2)^TA
651	//
652	ablascomplexsplitlength(ref a, m, ref s1, ref s2);
653	cmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
654	cmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2+s1, j2);
655	}
656	else
657	{
658
659	//
660	// Split A:
661	// (A1 A12)
662	// Xop(A) = Xop( )
663	// ( A2)
664	//
665	// Different variants depending on
666	// IsUpper/OpType combinations
667	//
668	ablascomplexsplitlength(ref a, n, ref s1, ref s2);
669	if( isupper & optype==0 )
670	{
671
672	//
673	// (A1 A12)-1
674	// XA^-1 = (X1 X2)( )
675	// ( A2)
676	//
677	cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
678	cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1+s1, 0, 1.0, ref x, i2, j2+s1);
679	cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
680	return;
681	}
682	if( isupper & optype!=0 )
683	{
684
685	//
686	// (A1' )-1
687	// XA^-1 = (X1 X2)( )
688	// (A12' A2')
689	//
690	cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
691	cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1, j1+s1, optype, 1.0, ref x, i2, j2);
692	cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
693	return;
694	}
695	if( !isupper & optype==0 )
696	{
697
698	//
699	// (A1 )-1
700	// XA^-1 = (X1 X2)( )
701	// (A21 A2)
702	//
703	cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
704	cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1+s1, j1, 0, 1.0, ref x, i2, j2);
705	cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
706	return;
707	}
708	if( !isupper & optype!=0 )
709	{
710
711	//
712	// (A1' A21')-1
713	// XA^-1 = (X1 X2)( )
714	// ( A2')
715	//
716	cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
717	cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1+s1, j1, optype, 1.0, ref x, i2, j2+s1);
718	cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
719	return;
720	}
721	}
722	}
723
724
725	/*************************************************************************
726	This subroutine calculates op(A^-1)*X where:
727	* X is MxN general matrix
728	* A is MxM upper/lower triangular/unitriangular matrix
729	* "op" may be identity transformation, transposition, conjugate transposition
730
731	Multiplication result replaces X.
732	Cache-oblivious algorithm is used.
733
734	INPUT PARAMETERS
735	N - matrix size, N>=0
736	M - matrix size, N>=0
737	A - matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1]
738	I1 - submatrix offset
739	J1 - submatrix offset
740	IsUpper - whether matrix is upper triangular
741	IsUnit - whether matrix is unitriangular
742	OpType - transformation type:
743	* 0 - no transformation
744	* 1 - transposition
745	* 2 - conjugate transposition
746	C - matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
747	I2 - submatrix offset
748	J2 - submatrix offset
749
750	-- ALGLIB routine --
751	15.12.2009
752	Bochkanov Sergey
753	*************************************************************************/
754	public static void cmatrixlefttrsm(int m,
755	int n,
756	ref AP.Complex[,] a,
757	int i1,
758	int j1,
759	bool isupper,
760	bool isunit,
761	int optype,
762	ref AP.Complex[,] x,
763	int i2,
764	int j2)
765	{
766	int s1 = 0;
767	int s2 = 0;
768	int bs = 0;
769
770	bs = ablascomplexblocksize(ref a);
771	if( m<=bs & n<=bs )
772	{
773	cmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
774	return;
775	}
776	if( n>=m )
777	{
778
779	//
780	// Split X: op(A)^-1X = op(A)^-1(X1 X2)
781	//
782	ablascomplexsplitlength(ref x, n, ref s1, ref s2);
783	cmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
784	cmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2+s1);
785	}
786	else
787	{
788
789	//
790	// Split A
791	//
792	ablascomplexsplitlength(ref a, m, ref s1, ref s2);
793	if( isupper & optype==0 )
794	{
795
796	//
797	// (A1 A12)-1 ( X1 )
798	// A^-1X = ( ) *( )
799	// ( A2) ( X2 )
800	//
801	cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
802	cmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1+s1, 0, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
803	cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
804	return;
805	}
806	if( isupper & optype!=0 )
807	{
808
809	//
810	// (A1' )-1 ( X1 )
811	// A^-1X = ( ) ( )
812	// (A12' A2') ( X2 )
813	//
814	cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
815	cmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1+s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
816	cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
817	return;
818	}
819	if( !isupper & optype==0 )
820	{
821
822	//
823	// (A1 )-1 ( X1 )
824	// A^-1X = ( ) ( )
825	// (A21 A2) ( X2 )
826	//
827	cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
828	cmatrixgemm(s2, n, s1, -1.0, ref a, i1+s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
829	cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
830	return;
831	}
832	if( !isupper & optype!=0 )
833	{
834
835	//
836	// (A1' A21')-1 ( X1 )
837	// A^-1X = ( ) ( )
838	// ( A2') ( X2 )
839	//
840	cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
841	cmatrixgemm(s1, n, s2, -1.0, ref a, i1+s1, j1, optype, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
842	cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
843	return;
844	}
845	}
846	}
847
848
849	/*************************************************************************
850	Same as CMatrixRightTRSM, but for real matrices
851
852	OpType may be only 0 or 1.
853
854	-- ALGLIB routine --
855	15.12.2009
856	Bochkanov Sergey
857	*************************************************************************/
858	public static void rmatrixrighttrsm(int m,
859	int n,
860	ref double[,] a,
861	int i1,
862	int j1,
863	bool isupper,
864	bool isunit,
865	int optype,
866	ref double[,] x,
867	int i2,
868	int j2)
869	{
870	int s1 = 0;
871	int s2 = 0;
872	int bs = 0;
873
874	bs = ablasblocksize(ref a);
875	if( m<=bs & n<=bs )
876	{
877	rmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
878	return;
879	}
880	if( m>=n )
881	{
882
883	//
884	// Split X: XA = (X1 X2)^TA
885	//
886	ablassplitlength(ref a, m, ref s1, ref s2);
887	rmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
888	rmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2+s1, j2);
889	}
890	else
891	{
892
893	//
894	// Split A:
895	// (A1 A12)
896	// Xop(A) = Xop( )
897	// ( A2)
898	//
899	// Different variants depending on
900	// IsUpper/OpType combinations
901	//
902	ablassplitlength(ref a, n, ref s1, ref s2);
903	if( isupper & optype==0 )
904	{
905
906	//
907	// (A1 A12)-1
908	// XA^-1 = (X1 X2)( )
909	// ( A2)
910	//
911	rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
912	rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1+s1, 0, 1.0, ref x, i2, j2+s1);
913	rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
914	return;
915	}
916	if( isupper & optype!=0 )
917	{
918
919	//
920	// (A1' )-1
921	// XA^-1 = (X1 X2)( )
922	// (A12' A2')
923	//
924	rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
925	rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1, j1+s1, optype, 1.0, ref x, i2, j2);
926	rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
927	return;
928	}
929	if( !isupper & optype==0 )
930	{
931
932	//
933	// (A1 )-1
934	// XA^-1 = (X1 X2)( )
935	// (A21 A2)
936	//
937	rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
938	rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1+s1, j1, 0, 1.0, ref x, i2, j2);
939	rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
940	return;
941	}
942	if( !isupper & optype!=0 )
943	{
944
945	//
946	// (A1' A21')-1
947	// XA^-1 = (X1 X2)( )
948	// ( A2')
949	//
950	rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
951	rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1+s1, j1, optype, 1.0, ref x, i2, j2+s1);
952	rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
953	return;
954	}
955	}
956	}
957
958
959	/*************************************************************************
960	Same as CMatrixLeftTRSM, but for real matrices
961
962	OpType may be only 0 or 1.
963
964	-- ALGLIB routine --
965	15.12.2009
966	Bochkanov Sergey
967	*************************************************************************/
968	public static void rmatrixlefttrsm(int m,
969	int n,
970	ref double[,] a,
971	int i1,
972	int j1,
973	bool isupper,
974	bool isunit,
975	int optype,
976	ref double[,] x,
977	int i2,
978	int j2)
979	{
980	int s1 = 0;
981	int s2 = 0;
982	int bs = 0;
983
984	bs = ablasblocksize(ref a);
985	if( m<=bs & n<=bs )
986	{
987	rmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
988	return;
989	}
990	if( n>=m )
991	{
992
993	//
994	// Split X: op(A)^-1X = op(A)^-1(X1 X2)
995	//
996	ablassplitlength(ref x, n, ref s1, ref s2);
997	rmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
998	rmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2+s1);
999	}
1000	else
1001	{
1002
1003	//
1004	// Split A
1005	//
1006	ablassplitlength(ref a, m, ref s1, ref s2);
1007	if( isupper & optype==0 )
1008	{
1009
1010	//
1011	// (A1 A12)-1 ( X1 )
1012	// A^-1X = ( ) *( )
1013	// ( A2) ( X2 )
1014	//
1015	rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
1016	rmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1+s1, 0, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
1017	rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
1018	return;
1019	}
1020	if( isupper & optype!=0 )
1021	{
1022
1023	//
1024	// (A1' )-1 ( X1 )
1025	// A^-1X = ( ) ( )
1026	// (A12' A2') ( X2 )
1027	//
1028	rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
1029	rmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1+s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
1030	rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
1031	return;
1032	}
1033	if( !isupper & optype==0 )
1034	{
1035
1036	//
1037	// (A1 )-1 ( X1 )
1038	// A^-1X = ( ) ( )
1039	// (A21 A2) ( X2 )
1040	//
1041	rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
1042	rmatrixgemm(s2, n, s1, -1.0, ref a, i1+s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
1043	rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
1044	return;
1045	}
1046	if( !isupper & optype!=0 )
1047	{
1048
1049	//
1050	// (A1' A21')-1 ( X1 )
1051	// A^-1X = ( ) ( )
1052	// ( A2') ( X2 )
1053	//
1054	rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
1055	rmatrixgemm(s1, n, s2, -1.0, ref a, i1+s1, j1, optype, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
1056	rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
1057	return;
1058	}
1059	}
1060	}
1061
1062
1063	/*************************************************************************
1064	This subroutine calculates C=alphaAA^H+betaC or C=alphaA^HA+betaC
1065	where:
1066	* C is NxN Hermitian matrix given by its upper/lower triangle
1067	* A is NxK matrix when A*A^H is calculated, KxN matrix otherwise
1068
1069	Additional info:
1070	* cache-oblivious algorithm is used.
1071	* multiplication result replaces C. If Beta=0, C elements are not used in
1072	calculations (not multiplied by zero - just not referenced)
1073	* if Alpha=0, A is not used (not multiplied by zero - just not referenced)
1074	* if both Beta and Alpha are zero, C is filled by zeros.
1075
1076	INPUT PARAMETERS
1077	N - matrix size, N>=0
1078	K - matrix size, K>=0
1079	Alpha - coefficient
1080	A - matrix
1081	IA - submatrix offset
1082	JA - submatrix offset
1083	OpTypeA - multiplication type:
1084	* 0 - A*A^H is calculated
1085	* 2 - A^H*A is calculated
1086	Beta - coefficient
1087	C - matrix
1088	IC - submatrix offset
1089	JC - submatrix offset
1090	IsUpper - whether C is upper triangular or lower triangular
1091
1092	-- ALGLIB routine --
1093	16.12.2009
1094	Bochkanov Sergey
1095	*************************************************************************/
1096	public static void cmatrixsyrk(int n,
1097	int k,
1098	double alpha,
1099	ref AP.Complex[,] a,
1100	int ia,
1101	int ja,
1102	int optypea,
1103	double beta,
1104	ref AP.Complex[,] c,
1105	int ic,
1106	int jc,
1107	bool isupper)
1108	{
1109	int s1 = 0;
1110	int s2 = 0;
1111	int bs = 0;
1112
1113	bs = ablascomplexblocksize(ref a);
1114	if( n<=bs & k<=bs )
1115	{
1116	cmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1117	return;
1118	}
1119	if( k>=n )
1120	{
1121
1122	//
1123	// Split K
1124	//
1125	ablascomplexsplitlength(ref a, k, ref s1, ref s2);
1126	if( optypea==0 )
1127	{
1128	cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1129	cmatrixsyrk(n, s2, alpha, ref a, ia, ja+s1, optypea, 1.0, ref c, ic, jc, isupper);
1130	}
1131	else
1132	{
1133	cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1134	cmatrixsyrk(n, s2, alpha, ref a, ia+s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
1135	}
1136	}
1137	else
1138	{
1139
1140	//
1141	// Split N
1142	//
1143	ablascomplexsplitlength(ref a, n, ref s1, ref s2);
1144	if( optypea==0 & isupper )
1145	{
1146	cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1147	cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia+s1, ja, 2, beta, ref c, ic, jc+s1);
1148	cmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1149	return;
1150	}
1151	if( optypea==0 & !isupper )
1152	{
1153	cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1154	cmatrixgemm(s2, s1, k, alpha, ref a, ia+s1, ja, 0, ref a, ia, ja, 2, beta, ref c, ic+s1, jc);
1155	cmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1156	return;
1157	}
1158	if( optypea!=0 & isupper )
1159	{
1160	cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1161	cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 2, ref a, ia, ja+s1, 0, beta, ref c, ic, jc+s1);
1162	cmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1163	return;
1164	}
1165	if( optypea!=0 & !isupper )
1166	{
1167	cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1168	cmatrixgemm(s2, s1, k, alpha, ref a, ia, ja+s1, 2, ref a, ia, ja, 0, beta, ref c, ic+s1, jc);
1169	cmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1170	return;
1171	}
1172	}
1173	}
1174
1175
1176	/*************************************************************************
1177	Same as CMatrixSYRK, but for real matrices
1178
1179	OpType may be only 0 or 1.
1180
1181	-- ALGLIB routine --
1182	16.12.2009
1183	Bochkanov Sergey
1184	*************************************************************************/
1185	public static void rmatrixsyrk(int n,
1186	int k,
1187	double alpha,
1188	ref double[,] a,
1189	int ia,
1190	int ja,
1191	int optypea,
1192	double beta,
1193	ref double[,] c,
1194	int ic,
1195	int jc,
1196	bool isupper)
1197	{
1198	int s1 = 0;
1199	int s2 = 0;
1200	int bs = 0;
1201
1202	bs = ablasblocksize(ref a);
1203	if( n<=bs & k<=bs )
1204	{
1205	rmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1206	return;
1207	}
1208	if( k>=n )
1209	{
1210
1211	//
1212	// Split K
1213	//
1214	ablassplitlength(ref a, k, ref s1, ref s2);
1215	if( optypea==0 )
1216	{
1217	rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1218	rmatrixsyrk(n, s2, alpha, ref a, ia, ja+s1, optypea, 1.0, ref c, ic, jc, isupper);
1219	}
1220	else
1221	{
1222	rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1223	rmatrixsyrk(n, s2, alpha, ref a, ia+s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
1224	}
1225	}
1226	else
1227	{
1228
1229	//
1230	// Split N
1231	//
1232	ablassplitlength(ref a, n, ref s1, ref s2);
1233	if( optypea==0 & isupper )
1234	{
1235	rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1236	rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia+s1, ja, 1, beta, ref c, ic, jc+s1);
1237	rmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1238	return;
1239	}
1240	if( optypea==0 & !isupper )
1241	{
1242	rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1243	rmatrixgemm(s2, s1, k, alpha, ref a, ia+s1, ja, 0, ref a, ia, ja, 1, beta, ref c, ic+s1, jc);
1244	rmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1245	return;
1246	}
1247	if( optypea!=0 & isupper )
1248	{
1249	rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1250	rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 1, ref a, ia, ja+s1, 0, beta, ref c, ic, jc+s1);
1251	rmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1252	return;
1253	}
1254	if( optypea!=0 & !isupper )
1255	{
1256	rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
1257	rmatrixgemm(s2, s1, k, alpha, ref a, ia, ja+s1, 1, ref a, ia, ja, 0, beta, ref c, ic+s1, jc);
1258	rmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
1259	return;
1260	}
1261	}
1262	}
1263
1264
1265	/*************************************************************************
1266	This subroutine calculates C = alphaop1(A)op2(B) +beta*C where:
1267	* C is MxN general matrix
1268	* op1(A) is MxK matrix
1269	* op2(B) is KxN matrix
1270	* "op" may be identity transformation, transposition, conjugate transposition
1271
1272	Additional info:
1273	* cache-oblivious algorithm is used.
1274	* multiplication result replaces C. If Beta=0, C elements are not used in
1275	calculations (not multiplied by zero - just not referenced)
1276	* if Alpha=0, A is not used (not multiplied by zero - just not referenced)
1277	* if both Beta and Alpha are zero, C is filled by zeros.
1278
1279	INPUT PARAMETERS
1280	N - matrix size, N>0
1281	M - matrix size, N>0
1282	K - matrix size, K>0
1283	Alpha - coefficient
1284	A - matrix
1285	IA - submatrix offset
1286	JA - submatrix offset
1287	OpTypeA - transformation type:
1288	* 0 - no transformation
1289	* 1 - transposition
1290	* 2 - conjugate transposition
1291	B - matrix
1292	IB - submatrix offset
1293	JB - submatrix offset
1294	OpTypeB - transformation type:
1295	* 0 - no transformation
1296	* 1 - transposition
1297	* 2 - conjugate transposition
1298	Beta - coefficient
1299	C - matrix
1300	IC - submatrix offset
1301	JC - submatrix offset
1302
1303	-- ALGLIB routine --
1304	16.12.2009
1305	Bochkanov Sergey
1306	*************************************************************************/
1307	public static void cmatrixgemm(int m,
1308	int n,
1309	int k,
1310	AP.Complex alpha,
1311	ref AP.Complex[,] a,
1312	int ia,
1313	int ja,
1314	int optypea,
1315	ref AP.Complex[,] b,
1316	int ib,
1317	int jb,
1318	int optypeb,
1319	AP.Complex beta,
1320	ref AP.Complex[,] c,
1321	int ic,
1322	int jc)
1323	{
1324	int s1 = 0;
1325	int s2 = 0;
1326	int bs = 0;
1327
1328	bs = ablascomplexblocksize(ref a);
1329	if( m<=bs & n<=bs & k<=bs )
1330	{
1331	cmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1332	return;
1333	}
1334	if( m>=n & m>=k )
1335	{
1336
1337	//
1338	// AB = (A1 A2)^TB
1339	//
1340	ablascomplexsplitlength(ref a, m, ref s1, ref s2);
1341	cmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1342	if( optypea==0 )
1343	{
1344	cmatrixgemm(s2, n, k, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
1345	}
1346	else
1347	{
1348	cmatrixgemm(s2, n, k, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
1349	}
1350	return;
1351	}
1352	if( n>=m & n>=k )
1353	{
1354
1355	//
1356	// AB = A(B1 B2)
1357	//
1358	ablascomplexsplitlength(ref a, n, ref s1, ref s2);
1359	if( optypeb==0 )
1360	{
1361	cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1362	cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb+s1, optypeb, beta, ref c, ic, jc+s1);
1363	}
1364	else
1365	{
1366	cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1367	cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib+s1, jb, optypeb, beta, ref c, ic, jc+s1);
1368	}
1369	return;
1370	}
1371	if( k>=m & k>=n )
1372	{
1373
1374	//
1375	// AB = (A1 A2)(B1 B2)^T
1376	//
1377	ablascomplexsplitlength(ref a, k, ref s1, ref s2);
1378	if( optypea==0 & optypeb==0 )
1379	{
1380	cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1381	cmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
1382	}
1383	if( optypea==0 & optypeb!=0 )
1384	{
1385	cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1386	cmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
1387	}
1388	if( optypea!=0 & optypeb==0 )
1389	{
1390	cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1391	cmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
1392	}
1393	if( optypea!=0 & optypeb!=0 )
1394	{
1395	cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1396	cmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
1397	}
1398	return;
1399	}
1400	}
1401
1402
1403	/*************************************************************************
1404	Same as CMatrixGEMM, but for real numbers.
1405	OpType may be only 0 or 1.
1406
1407	-- ALGLIB routine --
1408	16.12.2009
1409	Bochkanov Sergey
1410	*************************************************************************/
1411	public static void rmatrixgemm(int m,
1412	int n,
1413	int k,
1414	double alpha,
1415	ref double[,] a,
1416	int ia,
1417	int ja,
1418	int optypea,
1419	ref double[,] b,
1420	int ib,
1421	int jb,
1422	int optypeb,
1423	double beta,
1424	ref double[,] c,
1425	int ic,
1426	int jc)
1427	{
1428	int s1 = 0;
1429	int s2 = 0;
1430	int bs = 0;
1431
1432	bs = ablasblocksize(ref a);
1433	if( m<=bs & n<=bs & k<=bs )
1434	{
1435	rmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1436	return;
1437	}
1438	if( m>=n & m>=k )
1439	{
1440
1441	//
1442	// AB = (A1 A2)^TB
1443	//
1444	ablassplitlength(ref a, m, ref s1, ref s2);
1445	if( optypea==0 )
1446	{
1447	rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1448	rmatrixgemm(s2, n, k, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
1449	}
1450	else
1451	{
1452	rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1453	rmatrixgemm(s2, n, k, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
1454	}
1455	return;
1456	}
1457	if( n>=m & n>=k )
1458	{
1459
1460	//
1461	// AB = A(B1 B2)
1462	//
1463	ablassplitlength(ref a, n, ref s1, ref s2);
1464	if( optypeb==0 )
1465	{
1466	rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1467	rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb+s1, optypeb, beta, ref c, ic, jc+s1);
1468	}
1469	else
1470	{
1471	rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1472	rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib+s1, jb, optypeb, beta, ref c, ic, jc+s1);
1473	}
1474	return;
1475	}
1476	if( k>=m & k>=n )
1477	{
1478
1479	//
1480	// AB = (A1 A2)(B1 B2)^T
1481	//
1482	ablassplitlength(ref a, k, ref s1, ref s2);
1483	if( optypea==0 & optypeb==0 )
1484	{
1485	rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1486	rmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
1487	}
1488	if( optypea==0 & optypeb!=0 )
1489	{
1490	rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1491	rmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
1492	}
1493	if( optypea!=0 & optypeb==0 )
1494	{
1495	rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1496	rmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
1497	}
1498	if( optypea!=0 & optypeb!=0 )
1499	{
1500	rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
1501	rmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
1502	}
1503	return;
1504	}
1505	}
1506
1507
1508	/*************************************************************************
1509	Complex ABLASSplitLength
1510
1511	-- ALGLIB routine --
1512	15.12.2009
1513	Bochkanov Sergey
1514	*************************************************************************/
1515	private static void ablasinternalsplitlength(int n,
1516	int nb,
1517	ref int n1,
1518	ref int n2)
1519	{
1520	int r = 0;
1521
1522	if( n<=nb )
1523	{
1524
1525	//
1526	// Block size, no further splitting
1527	//
1528	n1 = n;
1529	n2 = 0;
1530	}
1531	else
1532	{
1533
1534	//
1535	// Greater than block size
1536	//
1537	if( n%nb!=0 )
1538	{
1539
1540	//
1541	// Split remainder
1542	//
1543	n2 = n%nb;
1544	n1 = n-n2;
1545	}
1546	else
1547	{
1548
1549	//
1550	// Split on block boundaries
1551	//
1552	n2 = n/2;
1553	n1 = n-n2;
1554	if( n1%nb==0 )
1555	{
1556	return;
1557	}
1558	r = nb-n1%nb;
1559	n1 = n1+r;
1560	n2 = n2-r;
1561	}
1562	}
1563	}
1564
1565
1566	/*************************************************************************
1567	Level 2 variant of CMatrixRightTRSM
1568	*************************************************************************/
1569	private static void cmatrixrighttrsm2(int m,
1570	int n,
1571	ref AP.Complex[,] a,
1572	int i1,
1573	int j1,
1574	bool isupper,
1575	bool isunit,
1576	int optype,
1577	ref AP.Complex[,] x,
1578	int i2,
1579	int j2)
1580	{
1581	int i = 0;
1582	int j = 0;
1583	AP.Complex vc = 0;
1584	AP.Complex vd = 0;
1585	int i_ = 0;
1586	int i1_ = 0;
1587
1588
1589	//
1590	// Special case
1591	//
1592	if( n*m==0 )
1593	{
1594	return;
1595	}
1596
1597	//
1598	// Try to call fast TRSM
1599	//
1600	if( ablasf.cmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
1601	{
1602	return;
1603	}
1604
1605	//
1606	// General case
1607	//
1608	if( isupper )
1609	{
1610
1611	//
1612	// Upper triangular matrix
1613	//
1614	if( optype==0 )
1615	{
1616
1617	//
1618	// X*A^(-1)
1619	//
1620	for(i=0; i<=m-1; i++)
1621	{
1622	for(j=0; j<=n-1; j++)
1623	{
1624	if( isunit )
1625	{
1626	vd = 1;
1627	}
1628	else
1629	{
1630	vd = a[i1+j,j1+j];
1631	}
1632	x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
1633	if( j<n-1 )
1634	{
1635	vc = x[i2+i,j2+j];
1636	i1_ = (j1+j+1) - (j2+j+1);
1637	for(i_=j2+j+1; i_<=j2+n-1;i_++)
1638	{
1639	x[i2+i,i_] = x[i2+i,i_] - vc*a[i1+j,i_+i1_];
1640	}
1641	}
1642	}
1643	}
1644	return;
1645	}
1646	if( optype==1 )
1647	{
1648
1649	//
1650	// X*A^(-T)
1651	//
1652	for(i=0; i<=m-1; i++)
1653	{
1654	for(j=n-1; j>=0; j--)
1655	{
1656	vc = 0;
1657	vd = 1;
1658	if( j<n-1 )
1659	{
1660	i1_ = (j1+j+1)-(j2+j+1);
1661	vc = 0.0;
1662	for(i_=j2+j+1; i_<=j2+n-1;i_++)
1663	{
1664	vc += x[i2+i,i_]*a[i1+j,i_+i1_];
1665	}
1666	}
1667	if( !isunit )
1668	{
1669	vd = a[i1+j,j1+j];
1670	}
1671	x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
1672	}
1673	}
1674	return;
1675	}
1676	if( optype==2 )
1677	{
1678
1679	//
1680	// X*A^(-H)
1681	//
1682	for(i=0; i<=m-1; i++)
1683	{
1684	for(j=n-1; j>=0; j--)
1685	{
1686	vc = 0;
1687	vd = 1;
1688	if( j<n-1 )
1689	{
1690	i1_ = (j1+j+1)-(j2+j+1);
1691	vc = 0.0;
1692	for(i_=j2+j+1; i_<=j2+n-1;i_++)
1693	{
1694	vc += x[i2+i,i_]*AP.Math.Conj(a[i1+j,i_+i1_]);
1695	}
1696	}
1697	if( !isunit )
1698	{
1699	vd = AP.Math.Conj(a[i1+j,j1+j]);
1700	}
1701	x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
1702	}
1703	}
1704	return;
1705	}
1706	}
1707	else
1708	{
1709
1710	//
1711	// Lower triangular matrix
1712	//
1713	if( optype==0 )
1714	{
1715
1716	//
1717	// X*A^(-1)
1718	//
1719	for(i=0; i<=m-1; i++)
1720	{
1721	for(j=n-1; j>=0; j--)
1722	{
1723	if( isunit )
1724	{
1725	vd = 1;
1726	}
1727	else
1728	{
1729	vd = a[i1+j,j1+j];
1730	}
1731	x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
1732	if( j>0 )
1733	{
1734	vc = x[i2+i,j2+j];
1735	i1_ = (j1) - (j2);
1736	for(i_=j2; i_<=j2+j-1;i_++)
1737	{
1738	x[i2+i,i_] = x[i2+i,i_] - vc*a[i1+j,i_+i1_];
1739	}
1740	}
1741	}
1742	}
1743	return;
1744	}
1745	if( optype==1 )
1746	{
1747
1748	//
1749	// X*A^(-T)
1750	//
1751	for(i=0; i<=m-1; i++)
1752	{
1753	for(j=0; j<=n-1; j++)
1754	{
1755	vc = 0;
1756	vd = 1;
1757	if( j>0 )
1758	{
1759	i1_ = (j1)-(j2);
1760	vc = 0.0;
1761	for(i_=j2; i_<=j2+j-1;i_++)
1762	{
1763	vc += x[i2+i,i_]*a[i1+j,i_+i1_];
1764	}
1765	}
1766	if( !isunit )
1767	{
1768	vd = a[i1+j,j1+j];
1769	}
1770	x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
1771	}
1772	}
1773	return;
1774	}
1775	if( optype==2 )
1776	{
1777
1778	//
1779	// X*A^(-H)
1780	//
1781	for(i=0; i<=m-1; i++)
1782	{
1783	for(j=0; j<=n-1; j++)
1784	{
1785	vc = 0;
1786	vd = 1;
1787	if( j>0 )
1788	{
1789	i1_ = (j1)-(j2);
1790	vc = 0.0;
1791	for(i_=j2; i_<=j2+j-1;i_++)
1792	{
1793	vc += x[i2+i,i_]*AP.Math.Conj(a[i1+j,i_+i1_]);
1794	}
1795	}
1796	if( !isunit )
1797	{
1798	vd = AP.Math.Conj(a[i1+j,j1+j]);
1799	}
1800	x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
1801	}
1802	}
1803	return;
1804	}
1805	}
1806	}
1807
1808
1809	/*************************************************************************
1810	Level-2 subroutine
1811	*************************************************************************/
1812	private static void cmatrixlefttrsm2(int m,
1813	int n,
1814	ref AP.Complex[,] a,
1815	int i1,
1816	int j1,
1817	bool isupper,
1818	bool isunit,
1819	int optype,
1820	ref AP.Complex[,] x,
1821	int i2,
1822	int j2)
1823	{
1824	int i = 0;
1825	int j = 0;
1826	AP.Complex vc = 0;
1827	AP.Complex vd = 0;
1828	int i_ = 0;
1829
1830
1831	//
1832	// Special case
1833	//
1834	if( n*m==0 )
1835	{
1836	return;
1837	}
1838
1839	//
1840	// Try to call fast TRSM
1841	//
1842	if( ablasf.cmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
1843	{
1844	return;
1845	}
1846
1847	//
1848	// General case
1849	//
1850	if( isupper )
1851	{
1852
1853	//
1854	// Upper triangular matrix
1855	//
1856	if( optype==0 )
1857	{
1858
1859	//
1860	// A^(-1)*X
1861	//
1862	for(i=m-1; i>=0; i--)
1863	{
1864	for(j=i+1; j<=m-1; j++)
1865	{
1866	vc = a[i1+i,j1+j];
1867	for(i_=j2; i_<=j2+n-1;i_++)
1868	{
1869	x[i2+i,i_] = x[i2+i,i_] - vc*x[i2+j,i_];
1870	}
1871	}
1872	if( !isunit )
1873	{
1874	vd = 1/a[i1+i,j1+i];
1875	for(i_=j2; i_<=j2+n-1;i_++)
1876	{
1877	x[i2+i,i_] = vd*x[i2+i,i_];
1878	}
1879	}
1880	}
1881	return;
1882	}
1883	if( optype==1 )
1884	{
1885
1886	//
1887	// A^(-T)*X
1888	//
1889	for(i=0; i<=m-1; i++)
1890	{
1891	if( isunit )
1892	{
1893	vd = 1;
1894	}
1895	else
1896	{
1897	vd = 1/a[i1+i,j1+i];
1898	}
1899	for(i_=j2; i_<=j2+n-1;i_++)
1900	{
1901	x[i2+i,i_] = vd*x[i2+i,i_];
1902	}
1903	for(j=i+1; j<=m-1; j++)
1904	{
1905	vc = a[i1+i,j1+j];
1906	for(i_=j2; i_<=j2+n-1;i_++)
1907	{
1908	x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
1909	}
1910	}
1911	}
1912	return;
1913	}
1914	if( optype==2 )
1915	{
1916
1917	//
1918	// A^(-H)*X
1919	//
1920	for(i=0; i<=m-1; i++)
1921	{
1922	if( isunit )
1923	{
1924	vd = 1;
1925	}
1926	else
1927	{
1928	vd = 1/AP.Math.Conj(a[i1+i,j1+i]);
1929	}
1930	for(i_=j2; i_<=j2+n-1;i_++)
1931	{
1932	x[i2+i,i_] = vd*x[i2+i,i_];
1933	}
1934	for(j=i+1; j<=m-1; j++)
1935	{
1936	vc = AP.Math.Conj(a[i1+i,j1+j]);
1937	for(i_=j2; i_<=j2+n-1;i_++)
1938	{
1939	x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
1940	}
1941	}
1942	}
1943	return;
1944	}
1945	}
1946	else
1947	{
1948
1949	//
1950	// Lower triangular matrix
1951	//
1952	if( optype==0 )
1953	{
1954
1955	//
1956	// A^(-1)*X
1957	//
1958	for(i=0; i<=m-1; i++)
1959	{
1960	for(j=0; j<=i-1; j++)
1961	{
1962	vc = a[i1+i,j1+j];
1963	for(i_=j2; i_<=j2+n-1;i_++)
1964	{
1965	x[i2+i,i_] = x[i2+i,i_] - vc*x[i2+j,i_];
1966	}
1967	}
1968	if( isunit )
1969	{
1970	vd = 1;
1971	}
1972	else
1973	{
1974	vd = 1/a[i1+j,j1+j];
1975	}
1976	for(i_=j2; i_<=j2+n-1;i_++)
1977	{
1978	x[i2+i,i_] = vd*x[i2+i,i_];
1979	}
1980	}
1981	return;
1982	}
1983	if( optype==1 )
1984	{
1985
1986	//
1987	// A^(-T)*X
1988	//
1989	for(i=m-1; i>=0; i--)
1990	{
1991	if( isunit )
1992	{
1993	vd = 1;
1994	}
1995	else
1996	{
1997	vd = 1/a[i1+i,j1+i];
1998	}
1999	for(i_=j2; i_<=j2+n-1;i_++)
2000	{
2001	x[i2+i,i_] = vd*x[i2+i,i_];
2002	}
2003	for(j=i-1; j>=0; j--)
2004	{
2005	vc = a[i1+i,j1+j];
2006	for(i_=j2; i_<=j2+n-1;i_++)
2007	{
2008	x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
2009	}
2010	}
2011	}
2012	return;
2013	}
2014	if( optype==2 )
2015	{
2016
2017	//
2018	// A^(-H)*X
2019	//
2020	for(i=m-1; i>=0; i--)
2021	{
2022	if( isunit )
2023	{
2024	vd = 1;
2025	}
2026	else
2027	{
2028	vd = 1/AP.Math.Conj(a[i1+i,j1+i]);
2029	}
2030	for(i_=j2; i_<=j2+n-1;i_++)
2031	{
2032	x[i2+i,i_] = vd*x[i2+i,i_];
2033	}
2034	for(j=i-1; j>=0; j--)
2035	{
2036	vc = AP.Math.Conj(a[i1+i,j1+j]);
2037	for(i_=j2; i_<=j2+n-1;i_++)
2038	{
2039	x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
2040	}
2041	}
2042	}
2043	return;
2044	}
2045	}
2046	}
2047
2048
2049	/*************************************************************************
2050	Level 2 subroutine
2051
2052	-- ALGLIB routine --
2053	15.12.2009
2054	Bochkanov Sergey
2055	*************************************************************************/
2056	private static void rmatrixrighttrsm2(int m,
2057	int n,
2058	ref double[,] a,
2059	int i1,
2060	int j1,
2061	bool isupper,
2062	bool isunit,
2063	int optype,
2064	ref double[,] x,
2065	int i2,
2066	int j2)
2067	{
2068	int i = 0;
2069	int j = 0;
2070	double vr = 0;
2071	double vd = 0;
2072	int i_ = 0;
2073	int i1_ = 0;
2074
2075
2076	//
2077	// Special case
2078	//
2079	if( n*m==0 )
2080	{
2081	return;
2082	}
2083
2084	//
2085	// Try to use "fast" code
2086	//
2087	if( ablasf.rmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
2088	{
2089	return;
2090	}
2091
2092	//
2093	// General case
2094	//
2095	if( isupper )
2096	{
2097
2098	//
2099	// Upper triangular matrix
2100	//
2101	if( optype==0 )
2102	{
2103
2104	//
2105	// X*A^(-1)
2106	//
2107	for(i=0; i<=m-1; i++)
2108	{
2109	for(j=0; j<=n-1; j++)
2110	{
2111	if( isunit )
2112	{
2113	vd = 1;
2114	}
2115	else
2116	{
2117	vd = a[i1+j,j1+j];
2118	}
2119	x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
2120	if( j<n-1 )
2121	{
2122	vr = x[i2+i,j2+j];
2123	i1_ = (j1+j+1) - (j2+j+1);
2124	for(i_=j2+j+1; i_<=j2+n-1;i_++)
2125	{
2126	x[i2+i,i_] = x[i2+i,i_] - vr*a[i1+j,i_+i1_];
2127	}
2128	}
2129	}
2130	}
2131	return;
2132	}
2133	if( optype==1 )
2134	{
2135
2136	//
2137	// X*A^(-T)
2138	//
2139	for(i=0; i<=m-1; i++)
2140	{
2141	for(j=n-1; j>=0; j--)
2142	{
2143	vr = 0;
2144	vd = 1;
2145	if( j<n-1 )
2146	{
2147	i1_ = (j1+j+1)-(j2+j+1);
2148	vr = 0.0;
2149	for(i_=j2+j+1; i_<=j2+n-1;i_++)
2150	{
2151	vr += x[i2+i,i_]*a[i1+j,i_+i1_];
2152	}
2153	}
2154	if( !isunit )
2155	{
2156	vd = a[i1+j,j1+j];
2157	}
2158	x[i2+i,j2+j] = (x[i2+i,j2+j]-vr)/vd;
2159	}
2160	}
2161	return;
2162	}
2163	}
2164	else
2165	{
2166
2167	//
2168	// Lower triangular matrix
2169	//
2170	if( optype==0 )
2171	{
2172
2173	//
2174	// X*A^(-1)
2175	//
2176	for(i=0; i<=m-1; i++)
2177	{
2178	for(j=n-1; j>=0; j--)
2179	{
2180	if( isunit )
2181	{
2182	vd = 1;
2183	}
2184	else
2185	{
2186	vd = a[i1+j,j1+j];
2187	}
2188	x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
2189	if( j>0 )
2190	{
2191	vr = x[i2+i,j2+j];
2192	i1_ = (j1) - (j2);
2193	for(i_=j2; i_<=j2+j-1;i_++)
2194	{
2195	x[i2+i,i_] = x[i2+i,i_] - vr*a[i1+j,i_+i1_];
2196	}
2197	}
2198	}
2199	}
2200	return;
2201	}
2202	if( optype==1 )
2203	{
2204
2205	//
2206	// X*A^(-T)
2207	//
2208	for(i=0; i<=m-1; i++)
2209	{
2210	for(j=0; j<=n-1; j++)
2211	{
2212	vr = 0;
2213	vd = 1;
2214	if( j>0 )
2215	{
2216	i1_ = (j1)-(j2);
2217	vr = 0.0;
2218	for(i_=j2; i_<=j2+j-1;i_++)
2219	{
2220	vr += x[i2+i,i_]*a[i1+j,i_+i1_];
2221	}
2222	}
2223	if( !isunit )
2224	{
2225	vd = a[i1+j,j1+j];
2226	}
2227	x[i2+i,j2+j] = (x[i2+i,j2+j]-vr)/vd;
2228	}
2229	}
2230	return;
2231	}
2232	}
2233	}
2234
2235
2236	/*************************************************************************
2237	Level 2 subroutine
2238	*************************************************************************/
2239	private static void rmatrixlefttrsm2(int m,
2240	int n,
2241	ref double[,] a,
2242	int i1,
2243	int j1,
2244	bool isupper,
2245	bool isunit,
2246	int optype,
2247	ref double[,] x,
2248	int i2,
2249	int j2)
2250	{
2251	int i = 0;
2252	int j = 0;
2253	double vr = 0;
2254	double vd = 0;
2255	int i_ = 0;
2256
2257
2258	//
2259	// Special case
2260	//
2261	if( n*m==0 )
2262	{
2263	return;
2264	}
2265
2266	//
2267	// Try fast code
2268	//
2269	if( ablasf.rmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
2270	{
2271	return;
2272	}
2273
2274	//
2275	// General case
2276	//
2277	if( isupper )
2278	{
2279
2280	//
2281	// Upper triangular matrix
2282	//
2283	if( optype==0 )
2284	{
2285
2286	//
2287	// A^(-1)*X
2288	//
2289	for(i=m-1; i>=0; i--)
2290	{
2291	for(j=i+1; j<=m-1; j++)
2292	{
2293	vr = a[i1+i,j1+j];
2294	for(i_=j2; i_<=j2+n-1;i_++)
2295	{
2296	x[i2+i,i_] = x[i2+i,i_] - vr*x[i2+j,i_];
2297	}
2298	}
2299	if( !isunit )
2300	{
2301	vd = 1/a[i1+i,j1+i];
2302	for(i_=j2; i_<=j2+n-1;i_++)
2303	{
2304	x[i2+i,i_] = vd*x[i2+i,i_];
2305	}
2306	}
2307	}
2308	return;
2309	}
2310	if( optype==1 )
2311	{
2312
2313	//
2314	// A^(-T)*X
2315	//
2316	for(i=0; i<=m-1; i++)
2317	{
2318	if( isunit )
2319	{
2320	vd = 1;
2321	}
2322	else
2323	{
2324	vd = 1/a[i1+i,j1+i];
2325	}
2326	for(i_=j2; i_<=j2+n-1;i_++)
2327	{
2328	x[i2+i,i_] = vd*x[i2+i,i_];
2329	}
2330	for(j=i+1; j<=m-1; j++)
2331	{
2332	vr = a[i1+i,j1+j];
2333	for(i_=j2; i_<=j2+n-1;i_++)
2334	{
2335	x[i2+j,i_] = x[i2+j,i_] - vr*x[i2+i,i_];
2336	}
2337	}
2338	}
2339	return;
2340	}
2341	}
2342	else
2343	{
2344
2345	//
2346	// Lower triangular matrix
2347	//
2348	if( optype==0 )
2349	{
2350
2351	//
2352	// A^(-1)*X
2353	//
2354	for(i=0; i<=m-1; i++)
2355	{
2356	for(j=0; j<=i-1; j++)
2357	{
2358	vr = a[i1+i,j1+j];
2359	for(i_=j2; i_<=j2+n-1;i_++)
2360	{
2361	x[i2+i,i_] = x[i2+i,i_] - vr*x[i2+j,i_];
2362	}
2363	}
2364	if( isunit )
2365	{
2366	vd = 1;
2367	}
2368	else
2369	{
2370	vd = 1/a[i1+j,j1+j];
2371	}
2372	for(i_=j2; i_<=j2+n-1;i_++)
2373	{
2374	x[i2+i,i_] = vd*x[i2+i,i_];
2375	}
2376	}
2377	return;
2378	}
2379	if( optype==1 )
2380	{
2381
2382	//
2383	// A^(-T)*X
2384	//
2385	for(i=m-1; i>=0; i--)
2386	{
2387	if( isunit )
2388	{
2389	vd = 1;
2390	}
2391	else
2392	{
2393	vd = 1/a[i1+i,j1+i];
2394	}
2395	for(i_=j2; i_<=j2+n-1;i_++)
2396	{
2397	x[i2+i,i_] = vd*x[i2+i,i_];
2398	}
2399	for(j=i-1; j>=0; j--)
2400	{
2401	vr = a[i1+i,j1+j];
2402	for(i_=j2; i_<=j2+n-1;i_++)
2403	{
2404	x[i2+j,i_] = x[i2+j,i_] - vr*x[i2+i,i_];
2405	}
2406	}
2407	}
2408	return;
2409	}
2410	}
2411	}
2412
2413
2414	/*************************************************************************
2415	Level 2 subroutine
2416	*************************************************************************/
2417	private static void cmatrixsyrk2(int n,
2418	int k,
2419	double alpha,
2420	ref AP.Complex[,] a,
2421	int ia,
2422	int ja,
2423	int optypea,
2424	double beta,
2425	ref AP.Complex[,] c,
2426	int ic,
2427	int jc,
2428	bool isupper)
2429	{
2430	int i = 0;
2431	int j = 0;
2432	int j1 = 0;
2433	int j2 = 0;
2434	AP.Complex v = 0;
2435	int i_ = 0;
2436	int i1_ = 0;
2437
2438
2439	//
2440	// Fast exit (nothing to be done)
2441	//
2442	if( ((double)(alpha)==(double)(0) \| k==0) & (double)(beta)==(double)(1) )
2443	{
2444	return;
2445	}
2446
2447	//
2448	// Try to call fast SYRK
2449	//
2450	if( ablasf.cmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper) )
2451	{
2452	return;
2453	}
2454
2455	//
2456	// SYRK
2457	//
2458	if( optypea==0 )
2459	{
2460
2461	//
2462	// C=alphaAA^H+beta*C
2463	//
2464	for(i=0; i<=n-1; i++)
2465	{
2466	if( isupper )
2467	{
2468	j1 = i;
2469	j2 = n-1;
2470	}
2471	else
2472	{
2473	j1 = 0;
2474	j2 = i;
2475	}
2476	for(j=j1; j<=j2; j++)
2477	{
2478	if( (double)(alpha)!=(double)(0) & k>0 )
2479	{
2480	v = 0.0;
2481	for(i_=ja; i_<=ja+k-1;i_++)
2482	{
2483	v += a[ia+i,i_]*AP.Math.Conj(a[ia+j,i_]);
2484	}
2485	}
2486	else
2487	{
2488	v = 0;
2489	}
2490	if( (double)(beta)==(double)(0) )
2491	{
2492	c[ic+i,jc+j] = alpha*v;
2493	}
2494	else
2495	{
2496	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
2497	}
2498	}
2499	}
2500	return;
2501	}
2502	else
2503	{
2504
2505	//
2506	// C=alphaA^HA+beta*C
2507	//
2508	for(i=0; i<=n-1; i++)
2509	{
2510	if( isupper )
2511	{
2512	j1 = i;
2513	j2 = n-1;
2514	}
2515	else
2516	{
2517	j1 = 0;
2518	j2 = i;
2519	}
2520	if( (double)(beta)==(double)(0) )
2521	{
2522	for(j=j1; j<=j2; j++)
2523	{
2524	c[ic+i,jc+j] = 0;
2525	}
2526	}
2527	else
2528	{
2529	for(i_=jc+j1; i_<=jc+j2;i_++)
2530	{
2531	c[ic+i,i_] = beta*c[ic+i,i_];
2532	}
2533	}
2534	}
2535	for(i=0; i<=k-1; i++)
2536	{
2537	for(j=0; j<=n-1; j++)
2538	{
2539	if( isupper )
2540	{
2541	j1 = j;
2542	j2 = n-1;
2543	}
2544	else
2545	{
2546	j1 = 0;
2547	j2 = j;
2548	}
2549	v = alpha*AP.Math.Conj(a[ia+i,ja+j]);
2550	i1_ = (ja+j1) - (jc+j1);
2551	for(i_=jc+j1; i_<=jc+j2;i_++)
2552	{
2553	c[ic+j,i_] = c[ic+j,i_] + v*a[ia+i,i_+i1_];
2554	}
2555	}
2556	}
2557	return;
2558	}
2559	}
2560
2561
2562	/*************************************************************************
2563	Level 2 subrotuine
2564	*************************************************************************/
2565	private static void rmatrixsyrk2(int n,
2566	int k,
2567	double alpha,
2568	ref double[,] a,
2569	int ia,
2570	int ja,
2571	int optypea,
2572	double beta,
2573	ref double[,] c,
2574	int ic,
2575	int jc,
2576	bool isupper)
2577	{
2578	int i = 0;
2579	int j = 0;
2580	int j1 = 0;
2581	int j2 = 0;
2582	double v = 0;
2583	int i_ = 0;
2584	int i1_ = 0;
2585
2586
2587	//
2588	// Fast exit (nothing to be done)
2589	//
2590	if( ((double)(alpha)==(double)(0) \| k==0) & (double)(beta)==(double)(1) )
2591	{
2592	return;
2593	}
2594
2595	//
2596	// Try to call fast SYRK
2597	//
2598	if( ablasf.rmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper) )
2599	{
2600	return;
2601	}
2602
2603	//
2604	// SYRK
2605	//
2606	if( optypea==0 )
2607	{
2608
2609	//
2610	// C=alphaAA^H+beta*C
2611	//
2612	for(i=0; i<=n-1; i++)
2613	{
2614	if( isupper )
2615	{
2616	j1 = i;
2617	j2 = n-1;
2618	}
2619	else
2620	{
2621	j1 = 0;
2622	j2 = i;
2623	}
2624	for(j=j1; j<=j2; j++)
2625	{
2626	if( (double)(alpha)!=(double)(0) & k>0 )
2627	{
2628	v = 0.0;
2629	for(i_=ja; i_<=ja+k-1;i_++)
2630	{
2631	v += a[ia+i,i_]*a[ia+j,i_];
2632	}
2633	}
2634	else
2635	{
2636	v = 0;
2637	}
2638	if( (double)(beta)==(double)(0) )
2639	{
2640	c[ic+i,jc+j] = alpha*v;
2641	}
2642	else
2643	{
2644	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
2645	}
2646	}
2647	}
2648	return;
2649	}
2650	else
2651	{
2652
2653	//
2654	// C=alphaA^HA+beta*C
2655	//
2656	for(i=0; i<=n-1; i++)
2657	{
2658	if( isupper )
2659	{
2660	j1 = i;
2661	j2 = n-1;
2662	}
2663	else
2664	{
2665	j1 = 0;
2666	j2 = i;
2667	}
2668	if( (double)(beta)==(double)(0) )
2669	{
2670	for(j=j1; j<=j2; j++)
2671	{
2672	c[ic+i,jc+j] = 0;
2673	}
2674	}
2675	else
2676	{
2677	for(i_=jc+j1; i_<=jc+j2;i_++)
2678	{
2679	c[ic+i,i_] = beta*c[ic+i,i_];
2680	}
2681	}
2682	}
2683	for(i=0; i<=k-1; i++)
2684	{
2685	for(j=0; j<=n-1; j++)
2686	{
2687	if( isupper )
2688	{
2689	j1 = j;
2690	j2 = n-1;
2691	}
2692	else
2693	{
2694	j1 = 0;
2695	j2 = j;
2696	}
2697	v = alpha*a[ia+i,ja+j];
2698	i1_ = (ja+j1) - (jc+j1);
2699	for(i_=jc+j1; i_<=jc+j2;i_++)
2700	{
2701	c[ic+j,i_] = c[ic+j,i_] + v*a[ia+i,i_+i1_];
2702	}
2703	}
2704	}
2705	return;
2706	}
2707	}
2708
2709
2710	/*************************************************************************
2711	GEMM kernel
2712
2713	-- ALGLIB routine --
2714	16.12.2009
2715	Bochkanov Sergey
2716	*************************************************************************/
2717	private static void cmatrixgemmk(int m,
2718	int n,
2719	int k,
2720	AP.Complex alpha,
2721	ref AP.Complex[,] a,
2722	int ia,
2723	int ja,
2724	int optypea,
2725	ref AP.Complex[,] b,
2726	int ib,
2727	int jb,
2728	int optypeb,
2729	AP.Complex beta,
2730	ref AP.Complex[,] c,
2731	int ic,
2732	int jc)
2733	{
2734	int i = 0;
2735	int j = 0;
2736	AP.Complex v = 0;
2737	int i_ = 0;
2738	int i1_ = 0;
2739
2740
2741	//
2742	// Special case
2743	//
2744	if( m*n==0 )
2745	{
2746	return;
2747	}
2748
2749	//
2750	// Try optimized code
2751	//
2752	if( ablasf.cmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc) )
2753	{
2754	return;
2755	}
2756
2757	//
2758	// Another special case
2759	//
2760	if( k==0 )
2761	{
2762	if( beta!=0 )
2763	{
2764	for(i=0; i<=m-1; i++)
2765	{
2766	for(j=0; j<=n-1; j++)
2767	{
2768	c[ic+i,jc+j] = beta*c[ic+i,jc+j];
2769	}
2770	}
2771	}
2772	else
2773	{
2774	for(i=0; i<=m-1; i++)
2775	{
2776	for(j=0; j<=n-1; j++)
2777	{
2778	c[ic+i,jc+j] = 0;
2779	}
2780	}
2781	}
2782	return;
2783	}
2784
2785	//
2786	// General case
2787	//
2788	if( optypea==0 & optypeb!=0 )
2789	{
2790
2791	//
2792	// A*B'
2793	//
2794	for(i=0; i<=m-1; i++)
2795	{
2796	for(j=0; j<=n-1; j++)
2797	{
2798	if( k==0 \| alpha==0 )
2799	{
2800	v = 0;
2801	}
2802	else
2803	{
2804	if( optypeb==1 )
2805	{
2806	i1_ = (jb)-(ja);
2807	v = 0.0;
2808	for(i_=ja; i_<=ja+k-1;i_++)
2809	{
2810	v += a[ia+i,i_]*b[ib+j,i_+i1_];
2811	}
2812	}
2813	else
2814	{
2815	i1_ = (jb)-(ja);
2816	v = 0.0;
2817	for(i_=ja; i_<=ja+k-1;i_++)
2818	{
2819	v += a[ia+i,i_]*AP.Math.Conj(b[ib+j,i_+i1_]);
2820	}
2821	}
2822	}
2823	if( beta==0 )
2824	{
2825	c[ic+i,jc+j] = alpha*v;
2826	}
2827	else
2828	{
2829	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
2830	}
2831	}
2832	}
2833	return;
2834	}
2835	if( optypea==0 & optypeb==0 )
2836	{
2837
2838	//
2839	// A*B
2840	//
2841	for(i=0; i<=m-1; i++)
2842	{
2843	if( beta!=0 )
2844	{
2845	for(i_=jc; i_<=jc+n-1;i_++)
2846	{
2847	c[ic+i,i_] = beta*c[ic+i,i_];
2848	}
2849	}
2850	else
2851	{
2852	for(j=0; j<=n-1; j++)
2853	{
2854	c[ic+i,jc+j] = 0;
2855	}
2856	}
2857	if( alpha!=0 )
2858	{
2859	for(j=0; j<=k-1; j++)
2860	{
2861	v = alpha*a[ia+i,ja+j];
2862	i1_ = (jb) - (jc);
2863	for(i_=jc; i_<=jc+n-1;i_++)
2864	{
2865	c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
2866	}
2867	}
2868	}
2869	}
2870	return;
2871	}
2872	if( optypea!=0 & optypeb!=0 )
2873	{
2874
2875	//
2876	// A'*B'
2877	//
2878	for(i=0; i<=m-1; i++)
2879	{
2880	for(j=0; j<=n-1; j++)
2881	{
2882	if( alpha==0 )
2883	{
2884	v = 0;
2885	}
2886	else
2887	{
2888	if( optypea==1 )
2889	{
2890	if( optypeb==1 )
2891	{
2892	i1_ = (jb)-(ia);
2893	v = 0.0;
2894	for(i_=ia; i_<=ia+k-1;i_++)
2895	{
2896	v += a[i_,ja+i]*b[ib+j,i_+i1_];
2897	}
2898	}
2899	else
2900	{
2901	i1_ = (jb)-(ia);
2902	v = 0.0;
2903	for(i_=ia; i_<=ia+k-1;i_++)
2904	{
2905	v += a[i_,ja+i]*AP.Math.Conj(b[ib+j,i_+i1_]);
2906	}
2907	}
2908	}
2909	else
2910	{
2911	if( optypeb==1 )
2912	{
2913	i1_ = (jb)-(ia);
2914	v = 0.0;
2915	for(i_=ia; i_<=ia+k-1;i_++)
2916	{
2917	v += AP.Math.Conj(a[i_,ja+i])*b[ib+j,i_+i1_];
2918	}
2919	}
2920	else
2921	{
2922	i1_ = (jb)-(ia);
2923	v = 0.0;
2924	for(i_=ia; i_<=ia+k-1;i_++)
2925	{
2926	v += AP.Math.Conj(a[i_,ja+i])*AP.Math.Conj(b[ib+j,i_+i1_]);
2927	}
2928	}
2929	}
2930	}
2931	if( beta==0 )
2932	{
2933	c[ic+i,jc+j] = alpha*v;
2934	}
2935	else
2936	{
2937	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
2938	}
2939	}
2940	}
2941	return;
2942	}
2943	if( optypea!=0 & optypeb==0 )
2944	{
2945
2946	//
2947	// A'*B
2948	//
2949	if( beta==0 )
2950	{
2951	for(i=0; i<=m-1; i++)
2952	{
2953	for(j=0; j<=n-1; j++)
2954	{
2955	c[ic+i,jc+j] = 0;
2956	}
2957	}
2958	}
2959	else
2960	{
2961	for(i=0; i<=m-1; i++)
2962	{
2963	for(i_=jc; i_<=jc+n-1;i_++)
2964	{
2965	c[ic+i,i_] = beta*c[ic+i,i_];
2966	}
2967	}
2968	}
2969	if( alpha!=0 )
2970	{
2971	for(j=0; j<=k-1; j++)
2972	{
2973	for(i=0; i<=m-1; i++)
2974	{
2975	if( optypea==1 )
2976	{
2977	v = alpha*a[ia+j,ja+i];
2978	}
2979	else
2980	{
2981	v = alpha*AP.Math.Conj(a[ia+j,ja+i]);
2982	}
2983	i1_ = (jb) - (jc);
2984	for(i_=jc; i_<=jc+n-1;i_++)
2985	{
2986	c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
2987	}
2988	}
2989	}
2990	}
2991	return;
2992	}
2993	}
2994
2995
2996	/*************************************************************************
2997	GEMM kernel
2998
2999	-- ALGLIB routine --
3000	16.12.2009
3001	Bochkanov Sergey
3002	*************************************************************************/
3003	private static void rmatrixgemmk(int m,
3004	int n,
3005	int k,
3006	double alpha,
3007	ref double[,] a,
3008	int ia,
3009	int ja,
3010	int optypea,
3011	ref double[,] b,
3012	int ib,
3013	int jb,
3014	int optypeb,
3015	double beta,
3016	ref double[,] c,
3017	int ic,
3018	int jc)
3019	{
3020	int i = 0;
3021	int j = 0;
3022	double v = 0;
3023	int i_ = 0;
3024	int i1_ = 0;
3025
3026
3027	//
3028	// if matrix size is zero
3029	//
3030	if( m*n==0 )
3031	{
3032	return;
3033	}
3034
3035	//
3036	// Try optimized code
3037	//
3038	if( ablasf.rmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc) )
3039	{
3040	return;
3041	}
3042
3043	//
3044	// if K=0, then C=Beta*C
3045	//
3046	if( k==0 )
3047	{
3048	if( (double)(beta)!=(double)(1) )
3049	{
3050	if( (double)(beta)!=(double)(0) )
3051	{
3052	for(i=0; i<=m-1; i++)
3053	{
3054	for(j=0; j<=n-1; j++)
3055	{
3056	c[ic+i,jc+j] = beta*c[ic+i,jc+j];
3057	}
3058	}
3059	}
3060	else
3061	{
3062	for(i=0; i<=m-1; i++)
3063	{
3064	for(j=0; j<=n-1; j++)
3065	{
3066	c[ic+i,jc+j] = 0;
3067	}
3068	}
3069	}
3070	}
3071	return;
3072	}
3073
3074	//
3075	// General case
3076	//
3077	if( optypea==0 & optypeb!=0 )
3078	{
3079
3080	//
3081	// A*B'
3082	//
3083	for(i=0; i<=m-1; i++)
3084	{
3085	for(j=0; j<=n-1; j++)
3086	{
3087	if( k==0 \| (double)(alpha)==(double)(0) )
3088	{
3089	v = 0;
3090	}
3091	else
3092	{
3093	i1_ = (jb)-(ja);
3094	v = 0.0;
3095	for(i_=ja; i_<=ja+k-1;i_++)
3096	{
3097	v += a[ia+i,i_]*b[ib+j,i_+i1_];
3098	}
3099	}
3100	if( (double)(beta)==(double)(0) )
3101	{
3102	c[ic+i,jc+j] = alpha*v;
3103	}
3104	else
3105	{
3106	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
3107	}
3108	}
3109	}
3110	return;
3111	}
3112	if( optypea==0 & optypeb==0 )
3113	{
3114
3115	//
3116	// A*B
3117	//
3118	for(i=0; i<=m-1; i++)
3119	{
3120	if( (double)(beta)!=(double)(0) )
3121	{
3122	for(i_=jc; i_<=jc+n-1;i_++)
3123	{
3124	c[ic+i,i_] = beta*c[ic+i,i_];
3125	}
3126	}
3127	else
3128	{
3129	for(j=0; j<=n-1; j++)
3130	{
3131	c[ic+i,jc+j] = 0;
3132	}
3133	}
3134	if( (double)(alpha)!=(double)(0) )
3135	{
3136	for(j=0; j<=k-1; j++)
3137	{
3138	v = alpha*a[ia+i,ja+j];
3139	i1_ = (jb) - (jc);
3140	for(i_=jc; i_<=jc+n-1;i_++)
3141	{
3142	c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
3143	}
3144	}
3145	}
3146	}
3147	return;
3148	}
3149	if( optypea!=0 & optypeb!=0 )
3150	{
3151
3152	//
3153	// A'*B'
3154	//
3155	for(i=0; i<=m-1; i++)
3156	{
3157	for(j=0; j<=n-1; j++)
3158	{
3159	if( (double)(alpha)==(double)(0) )
3160	{
3161	v = 0;
3162	}
3163	else
3164	{
3165	i1_ = (jb)-(ia);
3166	v = 0.0;
3167	for(i_=ia; i_<=ia+k-1;i_++)
3168	{
3169	v += a[i_,ja+i]*b[ib+j,i_+i1_];
3170	}
3171	}
3172	if( (double)(beta)==(double)(0) )
3173	{
3174	c[ic+i,jc+j] = alpha*v;
3175	}
3176	else
3177	{
3178	c[ic+i,jc+j] = betac[ic+i,jc+j]+alphav;
3179	}
3180	}
3181	}
3182	return;
3183	}
3184	if( optypea!=0 & optypeb==0 )
3185	{
3186
3187	//
3188	// A'*B
3189	//
3190	if( (double)(beta)==(double)(0) )
3191	{
3192	for(i=0; i<=m-1; i++)
3193	{
3194	for(j=0; j<=n-1; j++)
3195	{
3196	c[ic+i,jc+j] = 0;
3197	}
3198	}
3199	}
3200	else
3201	{
3202	for(i=0; i<=m-1; i++)
3203	{
3204	for(i_=jc; i_<=jc+n-1;i_++)
3205	{
3206	c[ic+i,i_] = beta*c[ic+i,i_];
3207	}
3208	}
3209	}
3210	if( (double)(alpha)!=(double)(0) )
3211	{
3212	for(j=0; j<=k-1; j++)
3213	{
3214	for(i=0; i<=m-1; i++)
3215	{
3216	v = alpha*a[ia+j,ja+i];
3217	i1_ = (jb) - (jc);
3218	for(i_=jc; i_<=jc+n-1;i_++)
3219	{
3220	c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
3221	}
3222	}
3223	}
3224	}
3225	return;
3226	}
3227	}
3228	}
3229	}

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

Download in other formats:

Update cookies preferences