Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/ablas.cs @ 3932

Visit:

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

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

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

Download in other formats:

Update cookies preferences