Context Navigation

source: branches/HeuristicLab.Problems.GaussianProcessTuning/ILNumerics.2.14.4735.573/Functions/builtin/eig.cs @ 11194

Visit:

Last change on this file since 11194 was 9102, checked in by gkronber, 12 years ago
#1967: ILNumerics source for experimentation
File size: 116.7 KB

Line
1	///
2	/// This file is part of ILNumerics Community Edition.
3	///
4	/// ILNumerics Community Edition - high performance computing for applications.
5	/// Copyright (C) 2006 - 2012 Haymo Kutschbach, http://ilnumerics.net
6	///
7	/// ILNumerics Community Edition is free software: you can redistribute it and/or modify
8	/// it under the terms of the GNU General Public License version 3 as published by
9	/// the Free Software Foundation.
10	///
11	/// ILNumerics Community Edition is distributed in the hope that it will be useful,
12	/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13	/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	/// GNU General Public License for more details.
15	///
16	/// You should have received a copy of the GNU General Public License
17	/// along with ILNumerics Community Edition. See the file License.txt in the root
18	/// of your distribution package. If not, see <http://www.gnu.org/licenses/>.
19	///
20	/// In addition this software uses the following components and/or licenses:
21	///
22	/// =================================================================================
23	/// The Open Toolkit Library License
24	///
25	/// Copyright (c) 2006 - 2009 the Open Toolkit library.
26	///
27	/// Permission is hereby granted, free of charge, to any person obtaining a copy
28	/// of this software and associated documentation files (the "Software"), to deal
29	/// in the Software without restriction, including without limitation the rights to
30	/// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
31	/// the Software, and to permit persons to whom the Software is furnished to do
32	/// so, subject to the following conditions:
33	///
34	/// The above copyright notice and this permission notice shall be included in all
35	/// copies or substantial portions of the Software.
36	///
37	/// =================================================================================
38	///
39
40	using System;
41	using System.Collections.Generic;
42	using System.Text;
43	using ILNumerics.Storage;
44	using ILNumerics.Misc;
45	using ILNumerics.Exceptions;
46
47
48
49	namespace ILNumerics {
50	public partial class ILMath {
51
52
53
54	/// <summary>
55	/// Compute eigenvalues of general square matrix A
56	/// </summary>
57	/// <param name="A">Input matrix A. Size [n x n]</param>
58	/// <returns>Vector of eigenvalues of A. Size [n x 1]</returns>
59	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
60	/// <para>The vector returned will be of complex inner type, since no further constraints are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
61	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
62	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>)"/>
63	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
64	public static ILRetArray< complex > eig(ILInArray< double > A) {
65	using (ILScope.Enter(A)) {
66	ILArray< complex> V = null;
67	MatrixProperties props = MatrixProperties.None;
68	return eig(A, V, ref props, true);
69	}
70	}
71	/// <summary>
72	/// Compute eigenvalues and eigenvectors of general square matrix A
73	/// </summary>
74	/// <param name="A">Input matrix A. Size [n x n]</param>
75	/// <param name="V">Output matrix, eigenvectors EV of size [n x n]. May be null on input. If not null, content of V will be destroyed.</param>
76	/// <returns>Diagonal matrix with eigenvalues of A. Size [n x n]</returns>
77	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
78	/// <para>The matrices returned will be of complex inner type, since no further constrains are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
79	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
80	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
81	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
82	public static ILRetArray< complex > eig(ILInArray< double > A, ILOutArray< complex > V) {
83	MatrixProperties props = MatrixProperties.None;
84	return eig(A, V, ref props, true);
85	}
86	/// <summary>
87	/// Find eigenvalues and eigenvectors
88	/// </summary>
89	/// <param name="A">Input: square matrix, size [n x n]</param>
90	/// <param name="V">Output (optional): eigenvectors</param>
91	/// <param name="propsA">Matrix properties, on input - if specified,
92	/// will be used to choose the proper method of solution. On exit will be
93	/// filled according to the properties of A.</param>
94	/// <param name="balance">true: permute A in order to increase the
95	/// numerical stability, false: do not permute A.</param>
96	/// <returns>eigenvalues as vector (if V is null) or as diagonoal
97	/// matrix (if V was requested, i.e. not null).</returns>
98	/// <remarks><para>The eigenvalues of A are found by use of the
99	/// Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
100	/// <para>The arrays returned will be of complex inner type,
101	/// since no further constraints are set on the structure of
102	/// A (it may be nonsymmetric). Use
103	/// <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/>
104	/// or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/>
105	/// functions for computing the real eigenvalues of symmetric
106	/// matrices explicitly.</para>
107	/// <para>Depending on the parameter <paramref name="balance"/>,
108	/// A will be balanced first. This includes permutations and
109	/// scaling of A in order to improve the conditioning of the
110	/// eigenvalues.</para></remarks>
111	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
112	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
113	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if a
114	/// is not square</exception>
115	public static ILRetArray< complex > eig(ILInArray< double > A, ILOutArray< complex > V, ref MatrixProperties propsA, bool balance) {
116	using (ILScope.Enter(A)) {
117	if (A.IsEmpty) {
118	if (!object.Equals(V,null))
119	V.a = empty<complex>(A.Size);
120	return empty<complex>(A.Size);
121	}
122	ILArray< complex > ret = empty< complex >(ILSize.Empty00);
123	int n = A.Size[0];
124	bool createVR = (object.Equals(V,null))? false:true;
125	if (n != A.Size[1])
126	throw new ILArgumentException("eig: matrix A must be square");
127	propsA \|= MatrixProperties.Square;
128	if (((propsA & MatrixProperties.Hermitian) != 0 \|\| ILMath.ishermitian(A.C))) {
129	propsA \|= MatrixProperties.Hermitian;
130	ILArray< double > Vd = null;
131	if (createVR)
132	Vd = returnType< double >();
133	ILArray< double > tmpRet = eigSymm(A,Vd);
134	if (createVR)
135	V.a = ILMath.tocomplex (Vd);
136	ret = ILMath.tocomplex (tmpRet);
137	} else {
138	// nonsymmetric case
139	char bal = (balance)? 'B':'N', jobvr;
140	ILRetArray< double > tmpA = A.C;
141	double [] vr = null;
142	double [] wr = ILMemoryPool.Pool.New< double >(n);
143
144	double[] wi = ILMemoryPool.Pool.New<double>(n);
145	double [] scale = ILMemoryPool.Pool.New< double >(n);
146	double [] rconde = ILMemoryPool.Pool.New< double >(n);
147	double [] rcondv = ILMemoryPool.Pool.New< double >(n);
148	double abnrm = 0;
149	int ldvr, ilo = 0, ihi = 0, info = 0;
150	if (createVR) {
151	ldvr = n;
152	vr = ILMemoryPool.Pool.New< double >(n * n);
153	jobvr = 'V';
154	} else {
155	ldvr = 1;
156	vr = new double [1];
157	jobvr = 'N';
158	}
159	/!HC:HC?geevx/
160	Lapack.dgeevx(bal,'N',jobvr,'N',n,tmpA.GetArrayForWrite(),n,wr,wi,new double [1],1,vr,ldvr,ref ilo,ref ihi,scale,ref abnrm,rconde,rcondv,ref info);
161	ILMemoryPool.Pool.Free(rconde);
162	ILMemoryPool.Pool.Free(rcondv);
163	ILMemoryPool.Pool.Free(scale);
164	if (info != 0)
165	throw new ILArgumentException("eig: error in Lapack '?geevx': (" + info + ")");
166	// create eigenvalues
167	complex [] retArr = ILMemoryPool.Pool.New< complex >(n);
168	for (int i = 0; i < n; i++) {
169
170	retArr[i].real = wr[i]; retArr[i].imag = wi[i];
171	}
172	ret = new ILArray< complex > (retArr,n,1);
173	if (createVR) {
174	#region HCSortEVec
175	complex [] VArr = ILMemoryPool.Pool.New< complex >(n * n);
176	for (int c = 0; c < n; c++) {
177	if (wi[c] != 0 && wi[c+1] != 0 && c < n-1) {
178	ilo = n * c; ihi = ilo + n;
179	for (int r = 0; r < n; r++) {
180	VArr[ilo].real = vr[ilo];
181	VArr[ilo].imag = vr[ihi];
182	VArr[ihi].real = vr[ilo];
183	VArr[ihi].imag = -vr[ihi];
184	ilo++; ihi++;
185	}
186	c++;
187	} else {
188	ilo = n * c;
189	for (int r = 0; r < n; r++) {
190	VArr[ilo].real = vr[ilo];
191	ilo++;
192	}
193	}
194	}
195	V.a = array<complex> (VArr,n,n);
196	#endregion HYCALPER
197	if (createVR)
198	ret.a = ILMath.diag< complex>(ret);
199	}
200
201	ILMemoryPool.Pool.Free(wi);
202	ILMemoryPool.Pool.Free(vr);
203	ILMemoryPool.Pool.Free(wr);
204	}
205	return ret;
206	}
207	}
208
209	/// <summary>
210	/// Find all eigenvalues of symmetric (hermitian) matrix
211	/// </summary>
212	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
213	/// <returns>Array of size [n,1] with eigenvalues of A.</returns>
214	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
215	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
216	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
217	public static ILRetArray< double > eigSymm (ILInArray< double > A) {
218	using (ILScope.Enter(A)) {
219	if (A.IsEmpty) {
220	return empty< double>(A.Size);
221	}
222	int n = A.Size[0];
223	if (n != A.Size[1])
224	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
225	int m = 0,info = 0;
226	ILArray< double > w = new ILArray< double > (new double [n],1,n);
227	double [] z = new double [1]; ;
228	int [] isuppz = new int[2 * n];
229	double [] AcArr = A.C.GetArrayForWrite();
230	/!HC:lapack_???evr/ Lapack.dsyevr ('N','A','U',n,AcArr,n,0,0,0,0,0,ref m,w.GetArrayForWrite(),z,1,isuppz,ref info);
231	ILMemoryPool.Pool.Free(AcArr);
232	return /dummy/ (w);
233	}
234	}
235	/// <summary>
236	/// Find all eigenvalues and -vectors of symmetric (hermitian) matrix
237	/// </summary>
238	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
239	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors
240	/// will not be computed and V is not changed. In order to make the function return the vectors, V should be initiialized with ILMath.returnType before calling eigSymm.</param>
241	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
242	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
243	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
244	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
245	public static ILRetArray< double > eigSymm (ILInArray< double > A, ILOutArray< double > V) {
246	using (ILScope.Enter(A)) {
247	if (A.IsEmpty) {
248	if (!object.Equals(V,null))
249	V.a = empty<double>(A.Size);
250	return empty<double>(A.Size);
251	}
252	int n = A.Size[0];
253	if (n != A.Size[1])
254	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
255	int m = 0,ldz = 0,info = 0;
256	ILArray< double > w = new ILArray< double >(ILMemoryPool.Pool.New< double >(n),n,1);
257	double [] z;
258	char jobz;
259	if (object.Equals(V,null)) {
260	z = new double [1];
261	jobz = 'N';
262	ldz = 1;
263	} else {
264	z = ILMemoryPool.Pool.New< double >(n * n);
265	jobz = 'V';
266	ldz = n;
267	}
268	int [] isuppz = ILMemoryPool.Pool.New<int>( 2 * n);
269	double [] AcArr = A.C.GetArrayForWrite();
270	/!HC:lapack_???evr/ Lapack.dsyevr (jobz,'A','U',n,AcArr,n,1,n,0,0,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
271	ILMemoryPool.Pool.Free(AcArr);
272	ILMemoryPool.Pool.Free(isuppz);
273	if (info != 0)
274	throw new ILException("error returned from lapack: " + info);
275	if (jobz == 'V') {
276	System.Diagnostics.Debug.Assert(!object.Equals(V,null));
277	V.a = array< double > (z,n,n);
278	V.a = V[full,r(0,m-1)];
279	return ILMath.diag( /dummy/ (w));
280	} else {
281	ILMemoryPool.Pool.Free(z);
282	return /dummy/ (w);
283	}
284	}
285	}
286	/// <summary>
287	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
288	/// </summary>
289	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
290	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
291	/// <param name="rangeStart">Specify the lowest limit for the range of eigenvalues to be queried.</param>
292	/// <param name="rangeEnd">Specify the upper limit for the range of eigenvalues to be queried.</param>
293	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
294	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
295	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
296	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/></exception>
297	public static ILRetArray< double > eigSymm (ILInArray< double > A, ILOutArray< double > V, int rangeStart, int rangeEnd) {
298	using (ILScope.Enter(A)) {
299	if (A.IsEmpty) {
300	if (!object.Equals(V,null))
301	V.a = empty<double>(A.Size);
302	return empty<double>(A.Size);
303	}
304	int n = A.Size[0];
305	if (n != A.Size[1])
306	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
307	int m = 0,ldz = 0,info = 0;
308	if (rangeEnd < rangeStart \|\| rangeStart < 1)
309	throw new ILArgumentException("invalid range of eigenvalues requested");
310	ILArray< double > w = array< double > (
311	ILMemoryPool.Pool.New< double>(n),1,n);
312	double [] z;
313	char jobz;
314	if (object.Equals(V,null)) {
315	z = new double [1];
316	jobz = 'N';
317	ldz = 1;
318	} else {
319	z = ILMemoryPool.Pool.New<double>(n * n);
320	jobz = 'V';
321	ldz = n;
322	}
323	int [] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
324	double [] AcArr = A.C.GetArrayForWrite();
325	/!HC:lapack_???evr/ Lapack.dsyevr (jobz,'I','U',n,AcArr,n,0,0,rangeStart,rangeEnd,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
326	ILMemoryPool.Pool.Free(isuppz);
327	ILMemoryPool.Pool.Free(AcArr);
328
329	if (jobz == 'V') {
330	V.a = array< double >(z,n,n);
331	V.a = V[full,r(0,m-1)];
332	}
333	ILMemoryPool.Pool.Free(z);
334	return ILMath.diag( /dummy/ (w));
335	}
336	}
337	/// <summary>
338	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
339	/// </summary>
340	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
341	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
342	/// <param name="rangeStart">The eigenvalues will be returned by increasing size. This will determine the number of the first eigenvalue to be returned.</param>
343	/// <param name="rangeEnd">Determine the number of the last eigenvalue to be returned.</param>
344	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
345	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
346	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
347	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/> or if either one is <= 0.</exception>
348	public static ILRetArray< double> eigSymm( ILInArray< double> A, ILOutArray< double> V, double rangeStart, double rangeEnd ) {
349	using (ILScope.Enter(A)) {
350	if (A.IsEmpty) {
351	if (!object.Equals(V, null))
352	V.a = empty<double>(A.Size);
353	return empty<double>(A.Size);
354	}
355	int n = A.Size[0];
356	if (n != A.Size[1])
357	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian");
358	int m = 0, ldz = 0, info = 0;
359	if (rangeStart > rangeEnd)
360	throw new ILArgumentException("invalid range of eigenvalues requested");
361	ILArray< double> w = zeros<double>(1, n);
362
363	double[] z;
364	char jobz;
365	if (object.Equals(V, null)) {
366	z = new double[1];
367	jobz = 'N';
368	ldz = 1;
369	} else {
370	z = ILMemoryPool.Pool.New<double>(n * n);
371	jobz = 'V';
372	ldz = n;
373	}
374	int[] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
375
376
377	double[] AcArr = A.C.GetArrayForWrite();
378	/!HC:lapack_???evr/
379	Lapack.dsyevr(jobz, 'V', 'U', n, AcArr, n, rangeStart, rangeEnd, 0, 0, 0, ref m, w.GetArrayForWrite(), z, ldz, isuppz, ref info);
380	ILMemoryPool.Pool.Free(AcArr);
381	ILMemoryPool.Pool.Free(isuppz);
382
383	if (jobz == 'V') {
384	V.a = array< double>(z, n, n);
385	V.a = V[full, r(0, m - 1)];
386	}
387	ILMemoryPool.Pool.Free(z);
388	return diag( /dummy/ (w));
389	}
390	}
391
392	#region HYCALPER AUTO GENERATED CODE
393
394
395	/// <summary>
396	/// Compute eigenvalues of general square matrix A
397	/// </summary>
398	/// <param name="A">Input matrix A. Size [n x n]</param>
399	/// <returns>Vector of eigenvalues of A. Size [n x 1]</returns>
400	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
401	/// <para>The vector returned will be of complex inner type, since no further constraints are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
402	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
403	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>)"/>
404	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
405	public static ILRetArray< fcomplex > eig(ILInArray< float > A) {
406	using (ILScope.Enter(A)) {
407	ILArray< fcomplex> V = null;
408	MatrixProperties props = MatrixProperties.None;
409	return eig(A, V, ref props, true);
410	}
411	}
412	/// <summary>
413	/// Compute eigenvalues and eigenvectors of general square matrix A
414	/// </summary>
415	/// <param name="A">Input matrix A. Size [n x n]</param>
416	/// <param name="V">Output matrix, eigenvectors EV of size [n x n]. May be null on input. If not null, content of V will be destroyed.</param>
417	/// <returns>Diagonal matrix with eigenvalues of A. Size [n x n]</returns>
418	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
419	/// <para>The matrices returned will be of complex inner type, since no further constrains are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
420	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
421	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
422	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
423	public static ILRetArray< fcomplex > eig(ILInArray< float > A, ILOutArray< fcomplex > V) {
424	MatrixProperties props = MatrixProperties.None;
425	return eig(A, V, ref props, true);
426	}
427	/// <summary>
428	/// Find eigenvalues and eigenvectors
429	/// </summary>
430	/// <param name="A">Input: square matrix, size [n x n]</param>
431	/// <param name="V">Output (optional): eigenvectors</param>
432	/// <param name="propsA">Matrix properties, on input - if specified,
433	/// will be used to choose the proper method of solution. On exit will be
434	/// filled according to the properties of A.</param>
435	/// <param name="balance">true: permute A in order to increase the
436	/// numerical stability, false: do not permute A.</param>
437	/// <returns>eigenvalues as vector (if V is null) or as diagonoal
438	/// matrix (if V was requested, i.e. not null).</returns>
439	/// <remarks><para>The eigenvalues of A are found by use of the
440	/// Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
441	/// <para>The arrays returned will be of complex inner type,
442	/// since no further constraints are set on the structure of
443	/// A (it may be nonsymmetric). Use
444	/// <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/>
445	/// or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/>
446	/// functions for computing the real eigenvalues of symmetric
447	/// matrices explicitly.</para>
448	/// <para>Depending on the parameter <paramref name="balance"/>,
449	/// A will be balanced first. This includes permutations and
450	/// scaling of A in order to improve the conditioning of the
451	/// eigenvalues.</para></remarks>
452	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
453	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
454	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if a
455	/// is not square</exception>
456	public static ILRetArray< fcomplex > eig(ILInArray< float > A, ILOutArray< fcomplex > V, ref MatrixProperties propsA, bool balance) {
457	using (ILScope.Enter(A)) {
458	if (A.IsEmpty) {
459	if (!object.Equals(V,null))
460	V.a = empty<fcomplex>(A.Size);
461	return empty<fcomplex>(A.Size);
462	}
463	ILArray< fcomplex > ret = empty< fcomplex >(ILSize.Empty00);
464	int n = A.Size[0];
465	bool createVR = (object.Equals(V,null))? false:true;
466	if (n != A.Size[1])
467	throw new ILArgumentException("eig: matrix A must be square");
468	propsA \|= MatrixProperties.Square;
469	if (((propsA & MatrixProperties.Hermitian) != 0 \|\| ILMath.ishermitian(A.C))) {
470	propsA \|= MatrixProperties.Hermitian;
471	ILArray< float > Vd = null;
472	if (createVR)
473	Vd = returnType< float >();
474	ILArray< float > tmpRet = eigSymm(A,Vd);
475	if (createVR)
476	V.a = ILMath.tofcomplex (Vd);
477	ret = ILMath.tofcomplex (tmpRet);
478	} else {
479	// nonsymmetric case
480	char bal = (balance)? 'B':'N', jobvr;
481	ILRetArray< float > tmpA = A.C;
482	float [] vr = null;
483	float [] wr = ILMemoryPool.Pool.New< float >(n);
484	float[] wi = ILMemoryPool.Pool.New< float>(n);
485	float [] scale = ILMemoryPool.Pool.New< float >(n);
486	float [] rconde = ILMemoryPool.Pool.New< float >(n);
487	float [] rcondv = ILMemoryPool.Pool.New< float >(n);
488	float abnrm = 0;
489	int ldvr, ilo = 0, ihi = 0, info = 0;
490	if (createVR) {
491	ldvr = n;
492	vr = ILMemoryPool.Pool.New< float >(n * n);
493	jobvr = 'V';
494	} else {
495	ldvr = 1;
496	vr = new float [1];
497	jobvr = 'N';
498	}
499	Lapack.sgeevx(bal,'N',jobvr,'N',n,tmpA.GetArrayForWrite(),n,wr,wi,new float [1],1,vr,ldvr,ref ilo,ref ihi,scale,ref abnrm,rconde,rcondv,ref info);
500	ILMemoryPool.Pool.Free(rconde);
501	ILMemoryPool.Pool.Free(rcondv);
502	ILMemoryPool.Pool.Free(scale);
503	if (info != 0)
504	throw new ILArgumentException("eig: error in Lapack '?geevx': (" + info + ")");
505	// create eigenvalues
506	fcomplex [] retArr = ILMemoryPool.Pool.New< fcomplex >(n);
507	for (int i = 0; i < n; i++) {
508	retArr[i].real = wr[i]; retArr[i].imag = wi[i];
509	}
510	ret = new ILArray< fcomplex > (retArr,n,1);
511	if (createVR) {
512	fcomplex [] VArr = ILMemoryPool.Pool.New< fcomplex >(n * n);
513	for (int c = 0; c < n; c++) {
514	if (wi[c] != 0 && wi[c+1] != 0 && c < n-1) {
515	ilo = n * c; ihi = ilo + n;
516	for (int r = 0; r < n; r++) {
517	VArr[ilo].real = vr[ilo];
518	VArr[ilo].imag = vr[ihi];
519	VArr[ihi].real = vr[ilo];
520	VArr[ihi].imag = -vr[ihi];
521	ilo++; ihi++;
522	}
523	c++;
524	} else {
525	ilo = n * c;
526	for (int r = 0; r < n; r++) {
527	VArr[ilo].real = vr[ilo];
528	ilo++;
529	}
530	}
531	}
532	V.a = array<fcomplex> (VArr,n,n);
533	if (createVR)
534	ret.a = ILMath.diag< fcomplex>(ret);
535	}
536	ILMemoryPool.Pool.Free(wi);
537	ILMemoryPool.Pool.Free(vr);
538	ILMemoryPool.Pool.Free(wr);
539	}
540	return ret;
541	}
542	}
543
544	/// <summary>
545	/// Find all eigenvalues of symmetric (hermitian) matrix
546	/// </summary>
547	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
548	/// <returns>Array of size [n,1] with eigenvalues of A.</returns>
549	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
550	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
551	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
552	public static ILRetArray< float > eigSymm (ILInArray< float > A) {
553	using (ILScope.Enter(A)) {
554	if (A.IsEmpty) {
555	return empty< float>(A.Size);
556	}
557	int n = A.Size[0];
558	if (n != A.Size[1])
559	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
560	int m = 0,info = 0;
561	ILArray< float > w = new ILArray< float > (new float [n],1,n);
562	float [] z = new float [1]; ;
563	int [] isuppz = new int[2 * n];
564	float [] AcArr = A.C.GetArrayForWrite();
565	Lapack.ssyevr ('N','A','U',n,AcArr,n,0,0,0,0,0,ref m,w.GetArrayForWrite(),z,1,isuppz,ref info);
566	ILMemoryPool.Pool.Free(AcArr);
567	return (w);
568	}
569	}
570	/// <summary>
571	/// Find all eigenvalues and -vectors of symmetric (hermitian) matrix
572	/// </summary>
573	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
574	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors
575	/// will not be computed and V is not changed. In order to make the function return the vectors, V should be initiialized with ILMath.returnType before calling eigSymm.</param>
576	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
577	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
578	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
579	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
580	public static ILRetArray< float > eigSymm (ILInArray< float > A, ILOutArray< float > V) {
581	using (ILScope.Enter(A)) {
582	if (A.IsEmpty) {
583	if (!object.Equals(V,null))
584	V.a = empty<float>(A.Size);
585	return empty<float>(A.Size);
586	}
587	int n = A.Size[0];
588	if (n != A.Size[1])
589	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
590	int m = 0,ldz = 0,info = 0;
591	ILArray< float > w = new ILArray< float >(ILMemoryPool.Pool.New< float >(n),n,1);
592	float [] z;
593	char jobz;
594	if (object.Equals(V,null)) {
595	z = new float [1];
596	jobz = 'N';
597	ldz = 1;
598	} else {
599	z = ILMemoryPool.Pool.New< float >(n * n);
600	jobz = 'V';
601	ldz = n;
602	}
603	int [] isuppz = ILMemoryPool.Pool.New<int>( 2 * n);
604	float [] AcArr = A.C.GetArrayForWrite();
605	Lapack.ssyevr (jobz,'A','U',n,AcArr,n,1,n,0,0,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
606	ILMemoryPool.Pool.Free(AcArr);
607	ILMemoryPool.Pool.Free(isuppz);
608	if (info != 0)
609	throw new ILException("error returned from lapack: " + info);
610	if (jobz == 'V') {
611	System.Diagnostics.Debug.Assert(!object.Equals(V,null));
612	V.a = array< float > (z,n,n);
613	V.a = V[full,r(0,m-1)];
614	return ILMath.diag( (w));
615	} else {
616	ILMemoryPool.Pool.Free(z);
617	return (w);
618	}
619	}
620	}
621	/// <summary>
622	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
623	/// </summary>
624	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
625	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
626	/// <param name="rangeStart">Specify the lowest limit for the range of eigenvalues to be queried.</param>
627	/// <param name="rangeEnd">Specify the upper limit for the range of eigenvalues to be queried.</param>
628	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
629	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
630	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
631	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/></exception>
632	public static ILRetArray< float > eigSymm (ILInArray< float > A, ILOutArray< float > V, int rangeStart, int rangeEnd) {
633	using (ILScope.Enter(A)) {
634	if (A.IsEmpty) {
635	if (!object.Equals(V,null))
636	V.a = empty<float>(A.Size);
637	return empty<float>(A.Size);
638	}
639	int n = A.Size[0];
640	if (n != A.Size[1])
641	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
642	int m = 0,ldz = 0,info = 0;
643	if (rangeEnd < rangeStart \|\| rangeStart < 1)
644	throw new ILArgumentException("invalid range of eigenvalues requested");
645	ILArray< float > w = array< float > (
646	ILMemoryPool.Pool.New< float>(n),1,n);
647	float [] z;
648	char jobz;
649	if (object.Equals(V,null)) {
650	z = new float [1];
651	jobz = 'N';
652	ldz = 1;
653	} else {
654	z = ILMemoryPool.Pool.New<float>(n * n);
655	jobz = 'V';
656	ldz = n;
657	}
658	int [] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
659	float [] AcArr = A.C.GetArrayForWrite();
660	Lapack.ssyevr (jobz,'I','U',n,AcArr,n,0,0,rangeStart,rangeEnd,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
661	ILMemoryPool.Pool.Free(isuppz);
662	ILMemoryPool.Pool.Free(AcArr);
663
664	if (jobz == 'V') {
665	V.a = array< float >(z,n,n);
666	V.a = V[full,r(0,m-1)];
667	}
668	ILMemoryPool.Pool.Free(z);
669	return ILMath.diag( (w));
670	}
671	}
672	/// <summary>
673	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
674	/// </summary>
675	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
676	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
677	/// <param name="rangeStart">The eigenvalues will be returned by increasing size. This will determine the number of the first eigenvalue to be returned.</param>
678	/// <param name="rangeEnd">Determine the number of the last eigenvalue to be returned.</param>
679	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
680	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
681	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
682	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/> or if either one is <= 0.</exception>
683	public static ILRetArray< float> eigSymm( ILInArray< float> A, ILOutArray< float> V, float rangeStart, float rangeEnd ) {
684	using (ILScope.Enter(A)) {
685	if (A.IsEmpty) {
686	if (!object.Equals(V, null))
687	V.a = empty<float>(A.Size);
688	return empty<float>(A.Size);
689	}
690	int n = A.Size[0];
691	if (n != A.Size[1])
692	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian");
693	int m = 0, ldz = 0, info = 0;
694	if (rangeStart > rangeEnd)
695	throw new ILArgumentException("invalid range of eigenvalues requested");
696	ILArray< float> w = zeros<float>(1, n);
697
698	float[] z;
699	char jobz;
700	if (object.Equals(V, null)) {
701	z = new float[1];
702	jobz = 'N';
703	ldz = 1;
704	} else {
705	z = ILMemoryPool.Pool.New<float>(n * n);
706	jobz = 'V';
707	ldz = n;
708	}
709	int[] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
710
711
712	float[] AcArr = A.C.GetArrayForWrite();
713
714	Lapack.ssyevr(jobz, 'V', 'U', n, AcArr, n, rangeStart, rangeEnd, 0, 0, 0, ref m, w.GetArrayForWrite(), z, ldz, isuppz, ref info);
715	ILMemoryPool.Pool.Free(AcArr);
716	ILMemoryPool.Pool.Free(isuppz);
717
718	if (jobz == 'V') {
719	V.a = array< float>(z, n, n);
720	V.a = V[full, r(0, m - 1)];
721	}
722	ILMemoryPool.Pool.Free(z);
723	return diag( (w));
724	}
725	}
726
727	/// <summary>
728	/// Compute eigenvalues of general square matrix A
729	/// </summary>
730	/// <param name="A">Input matrix A. Size [n x n]</param>
731	/// <returns>Vector of eigenvalues of A. Size [n x 1]</returns>
732	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
733	/// <para>The vector returned will be of complex inner type, since no further constraints are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
734	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
735	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>)"/>
736	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
737	public static ILRetArray< fcomplex > eig(ILInArray< fcomplex > A) {
738	using (ILScope.Enter(A)) {
739	ILArray< fcomplex> V = null;
740	MatrixProperties props = MatrixProperties.None;
741	return eig(A, V, ref props, true);
742	}
743	}
744	/// <summary>
745	/// Compute eigenvalues and eigenvectors of general square matrix A
746	/// </summary>
747	/// <param name="A">Input matrix A. Size [n x n]</param>
748	/// <param name="V">Output matrix, eigenvectors EV of size [n x n]. May be null on input. If not null, content of V will be destroyed.</param>
749	/// <returns>Diagonal matrix with eigenvalues of A. Size [n x n]</returns>
750	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
751	/// <para>The matrices returned will be of complex inner type, since no further constrains are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
752	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
753	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
754	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
755	public static ILRetArray< fcomplex > eig(ILInArray< fcomplex > A, ILOutArray< fcomplex > V) {
756	MatrixProperties props = MatrixProperties.None;
757	return eig(A, V, ref props, true);
758	}
759	/// <summary>
760	/// Find eigenvalues and eigenvectors
761	/// </summary>
762	/// <param name="A">Input: square matrix, size [n x n]</param>
763	/// <param name="V">Output (optional): eigenvectors</param>
764	/// <param name="propsA">Matrix properties, on input - if specified,
765	/// will be used to choose the proper method of solution. On exit will be
766	/// filled according to the properties of A.</param>
767	/// <param name="balance">true: permute A in order to increase the
768	/// numerical stability, false: do not permute A.</param>
769	/// <returns>eigenvalues as vector (if V is null) or as diagonoal
770	/// matrix (if V was requested, i.e. not null).</returns>
771	/// <remarks><para>The eigenvalues of A are found by use of the
772	/// Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
773	/// <para>The arrays returned will be of complex inner type,
774	/// since no further constraints are set on the structure of
775	/// A (it may be nonsymmetric). Use
776	/// <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/>
777	/// or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/>
778	/// functions for computing the real eigenvalues of symmetric
779	/// matrices explicitly.</para>
780	/// <para>Depending on the parameter <paramref name="balance"/>,
781	/// A will be balanced first. This includes permutations and
782	/// scaling of A in order to improve the conditioning of the
783	/// eigenvalues.</para></remarks>
784	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
785	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
786	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if a
787	/// is not square</exception>
788	public static ILRetArray< fcomplex > eig(ILInArray< fcomplex > A, ILOutArray< fcomplex > V, ref MatrixProperties propsA, bool balance) {
789	using (ILScope.Enter(A)) {
790	if (A.IsEmpty) {
791	if (!object.Equals(V,null))
792	V.a = empty<fcomplex>(A.Size);
793	return empty<fcomplex>(A.Size);
794	}
795	ILArray< fcomplex > ret = empty< fcomplex >(ILSize.Empty00);
796	int n = A.Size[0];
797	bool createVR = (object.Equals(V,null))? false:true;
798	if (n != A.Size[1])
799	throw new ILArgumentException("eig: matrix A must be square");
800	propsA \|= MatrixProperties.Square;
801	if (((propsA & MatrixProperties.Hermitian) != 0 \|\| ILMath.ishermitian(A.C))) {
802	propsA \|= MatrixProperties.Hermitian;
803	ILArray< fcomplex > Vd = null;
804	if (createVR)
805	Vd = returnType< fcomplex >();
806	ILArray< fcomplex > tmpRet = eigSymm(A,Vd);
807	if (createVR)
808	V.a = (Vd);
809	ret = (tmpRet);
810	} else {
811	// nonsymmetric case
812	char bal = (balance)? 'B':'N', jobvr;
813	ILRetArray< fcomplex > tmpA = A.C;
814	fcomplex [] vr = null;
815	fcomplex [] wr = ILMemoryPool.Pool.New< fcomplex >(n);
816
817	float [] scale = ILMemoryPool.Pool.New< float >(n);
818	float [] rconde = ILMemoryPool.Pool.New< float >(n);
819	float [] rcondv = ILMemoryPool.Pool.New< float >(n);
820	float abnrm = 0;
821	int ldvr, ilo = 0, ihi = 0, info = 0;
822	if (createVR) {
823	ldvr = n;
824	vr = ILMemoryPool.Pool.New< fcomplex >(n * n);
825	jobvr = 'V';
826	} else {
827	ldvr = 1;
828	vr = new fcomplex [1];
829	jobvr = 'N';
830	}
831	Lapack.cgeevx(bal,'N',jobvr,'N',n,tmpA.GetArrayForWrite(),n,wr, new fcomplex[1],1,vr,ldvr,ref ilo,ref ihi,scale,ref abnrm,rconde,rcondv,ref info);
832	ILMemoryPool.Pool.Free(rconde);
833	ILMemoryPool.Pool.Free(rcondv);
834	ILMemoryPool.Pool.Free(scale);
835	if (info != 0)
836	throw new ILArgumentException("eig: error in Lapack '?geevx': (" + info + ")");
837	// create eigenvalues
838	fcomplex [] retArr = ILMemoryPool.Pool.New< fcomplex >(n);
839	for (int i = 0; i < n; i++) {
840	retArr[i] = wr[i];
841	}
842	ret = new ILArray< fcomplex > (retArr,n,1);
843	if (createVR) {
844	V.a = array<fcomplex> (vr,n,n);
845	if (createVR)
846	ret.a = ILMath.diag< fcomplex>(ret);
847	}
848
849	ILMemoryPool.Pool.Free(vr);
850	ILMemoryPool.Pool.Free(wr);
851	}
852	return ret;
853	}
854	}
855
856	/// <summary>
857	/// Find all eigenvalues of symmetric (hermitian) matrix
858	/// </summary>
859	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
860	/// <returns>Array of size [n,1] with eigenvalues of A.</returns>
861	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
862	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
863	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
864	public static ILRetArray< fcomplex > eigSymm (ILInArray< fcomplex > A) {
865	using (ILScope.Enter(A)) {
866	if (A.IsEmpty) {
867	return empty< fcomplex>(A.Size);
868	}
869	int n = A.Size[0];
870	if (n != A.Size[1])
871	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
872	int m = 0,info = 0;
873	ILArray< float > w = new ILArray< float > (new float [n],1,n);
874	fcomplex [] z = new fcomplex [1]; ;
875	int [] isuppz = new int[2 * n];
876	fcomplex [] AcArr = A.C.GetArrayForWrite();
877	Lapack.cheevr ('N','A','U',n,AcArr,n,0,0,0,0,0,ref m,w.GetArrayForWrite(),z,1,isuppz,ref info);
878	ILMemoryPool.Pool.Free(AcArr);
879	return ILMath.tofcomplex(w);
880	}
881	}
882	/// <summary>
883	/// Find all eigenvalues and -vectors of symmetric (hermitian) matrix
884	/// </summary>
885	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
886	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors
887	/// will not be computed and V is not changed. In order to make the function return the vectors, V should be initiialized with ILMath.returnType before calling eigSymm.</param>
888	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
889	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
890	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
891	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
892	public static ILRetArray< fcomplex > eigSymm (ILInArray< fcomplex > A, ILOutArray< fcomplex > V) {
893	using (ILScope.Enter(A)) {
894	if (A.IsEmpty) {
895	if (!object.Equals(V,null))
896	V.a = empty<fcomplex>(A.Size);
897	return empty<fcomplex>(A.Size);
898	}
899	int n = A.Size[0];
900	if (n != A.Size[1])
901	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
902	int m = 0,ldz = 0,info = 0;
903	ILArray< float > w = new ILArray< float >(ILMemoryPool.Pool.New< float >(n),n,1);
904	fcomplex [] z;
905	char jobz;
906	if (object.Equals(V,null)) {
907	z = new fcomplex [1];
908	jobz = 'N';
909	ldz = 1;
910	} else {
911	z = ILMemoryPool.Pool.New< fcomplex >(n * n);
912	jobz = 'V';
913	ldz = n;
914	}
915	int [] isuppz = ILMemoryPool.Pool.New<int>( 2 * n);
916	fcomplex [] AcArr = A.C.GetArrayForWrite();
917	Lapack.cheevr (jobz,'A','U',n,AcArr,n,1,n,0,0,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
918	ILMemoryPool.Pool.Free(AcArr);
919	ILMemoryPool.Pool.Free(isuppz);
920	if (info != 0)
921	throw new ILException("error returned from lapack: " + info);
922	if (jobz == 'V') {
923	System.Diagnostics.Debug.Assert(!object.Equals(V,null));
924	V.a = array< fcomplex > (z,n,n);
925	V.a = V[full,r(0,m-1)];
926	return ILMath.diag( ILMath.tofcomplex(w));
927	} else {
928	ILMemoryPool.Pool.Free(z);
929	return ILMath.tofcomplex(w);
930	}
931	}
932	}
933	/// <summary>
934	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
935	/// </summary>
936	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
937	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
938	/// <param name="rangeStart">Specify the lowest limit for the range of eigenvalues to be queried.</param>
939	/// <param name="rangeEnd">Specify the upper limit for the range of eigenvalues to be queried.</param>
940	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
941	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
942	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
943	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/></exception>
944	public static ILRetArray< fcomplex > eigSymm (ILInArray< fcomplex > A, ILOutArray< fcomplex > V, int rangeStart, int rangeEnd) {
945	using (ILScope.Enter(A)) {
946	if (A.IsEmpty) {
947	if (!object.Equals(V,null))
948	V.a = empty<fcomplex>(A.Size);
949	return empty<fcomplex>(A.Size);
950	}
951	int n = A.Size[0];
952	if (n != A.Size[1])
953	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
954	int m = 0,ldz = 0,info = 0;
955	if (rangeEnd < rangeStart \|\| rangeStart < 1)
956	throw new ILArgumentException("invalid range of eigenvalues requested");
957	ILArray< float > w = array< float > (
958	ILMemoryPool.Pool.New< float>(n),1,n);
959	fcomplex [] z;
960	char jobz;
961	if (object.Equals(V,null)) {
962	z = new fcomplex [1];
963	jobz = 'N';
964	ldz = 1;
965	} else {
966	z = ILMemoryPool.Pool.New<fcomplex>(n * n);
967	jobz = 'V';
968	ldz = n;
969	}
970	int [] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
971	fcomplex [] AcArr = A.C.GetArrayForWrite();
972	Lapack.cheevr (jobz,'I','U',n,AcArr,n,0,0,rangeStart,rangeEnd,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
973	ILMemoryPool.Pool.Free(isuppz);
974	ILMemoryPool.Pool.Free(AcArr);
975
976	if (jobz == 'V') {
977	V.a = array< fcomplex >(z,n,n);
978	V.a = V[full,r(0,m-1)];
979	}
980	ILMemoryPool.Pool.Free(z);
981	return ILMath.diag( ILMath.tofcomplex(w));
982	}
983	}
984	/// <summary>
985	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
986	/// </summary>
987	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
988	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
989	/// <param name="rangeStart">The eigenvalues will be returned by increasing size. This will determine the number of the first eigenvalue to be returned.</param>
990	/// <param name="rangeEnd">Determine the number of the last eigenvalue to be returned.</param>
991	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
992	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
993	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
994	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/> or if either one is <= 0.</exception>
995	public static ILRetArray< fcomplex> eigSymm( ILInArray< fcomplex> A, ILOutArray< fcomplex> V, float rangeStart, float rangeEnd ) {
996	using (ILScope.Enter(A)) {
997	if (A.IsEmpty) {
998	if (!object.Equals(V, null))
999	V.a = empty<fcomplex>(A.Size);
1000	return empty<fcomplex>(A.Size);
1001	}
1002	int n = A.Size[0];
1003	if (n != A.Size[1])
1004	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian");
1005	int m = 0, ldz = 0, info = 0;
1006	if (rangeStart > rangeEnd)
1007	throw new ILArgumentException("invalid range of eigenvalues requested");
1008	ILArray< float> w = zeros<float>(1, n);
1009
1010	fcomplex[] z;
1011	char jobz;
1012	if (object.Equals(V, null)) {
1013	z = new fcomplex[1];
1014	jobz = 'N';
1015	ldz = 1;
1016	} else {
1017	z = ILMemoryPool.Pool.New<fcomplex>(n * n);
1018	jobz = 'V';
1019	ldz = n;
1020	}
1021	int[] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
1022
1023
1024	fcomplex[] AcArr = A.C.GetArrayForWrite();
1025
1026	Lapack.cheevr(jobz, 'V', 'U', n, AcArr, n, rangeStart, rangeEnd, 0, 0, 0, ref m, w.GetArrayForWrite(), z, ldz, isuppz, ref info);
1027	ILMemoryPool.Pool.Free(AcArr);
1028	ILMemoryPool.Pool.Free(isuppz);
1029
1030	if (jobz == 'V') {
1031	V.a = array< fcomplex>(z, n, n);
1032	V.a = V[full, r(0, m - 1)];
1033	}
1034	ILMemoryPool.Pool.Free(z);
1035	return diag( ILMath.tofcomplex(w));
1036	}
1037	}
1038
1039	/// <summary>
1040	/// Compute eigenvalues of general square matrix A
1041	/// </summary>
1042	/// <param name="A">Input matrix A. Size [n x n]</param>
1043	/// <returns>Vector of eigenvalues of A. Size [n x 1]</returns>
1044	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
1045	/// <para>The vector returned will be of complex inner type, since no further constraints are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
1046	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
1047	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>)"/>
1048	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
1049	public static ILRetArray< complex > eig(ILInArray< complex > A) {
1050	using (ILScope.Enter(A)) {
1051	ILArray< complex> V = null;
1052	MatrixProperties props = MatrixProperties.None;
1053	return eig(A, V, ref props, true);
1054	}
1055	}
1056	/// <summary>
1057	/// Compute eigenvalues and eigenvectors of general square matrix A
1058	/// </summary>
1059	/// <param name="A">Input matrix A. Size [n x n]</param>
1060	/// <param name="V">Output matrix, eigenvectors EV of size [n x n]. May be null on input. If not null, content of V will be destroyed.</param>
1061	/// <returns>Diagonal matrix with eigenvalues of A. Size [n x n]</returns>
1062	/// <remarks><para>The eigenvalues of A are found by use of the Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
1063	/// <para>The matrices returned will be of complex inner type, since no further constrains are set on the structure of A (it may be nonsymmetric). Use <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/> or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/> functions for computing the real eigenvalues of symmetric matrices explicitly.</para>
1064	/// <para>A will be balanced first. This includes permutations and scaling of A in order to improve the conditioning of the eigenvalues.</para></remarks>
1065	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
1066	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
1067	public static ILRetArray< complex > eig(ILInArray< complex > A, ILOutArray< complex > V) {
1068	MatrixProperties props = MatrixProperties.None;
1069	return eig(A, V, ref props, true);
1070	}
1071	/// <summary>
1072	/// Find eigenvalues and eigenvectors
1073	/// </summary>
1074	/// <param name="A">Input: square matrix, size [n x n]</param>
1075	/// <param name="V">Output (optional): eigenvectors</param>
1076	/// <param name="propsA">Matrix properties, on input - if specified,
1077	/// will be used to choose the proper method of solution. On exit will be
1078	/// filled according to the properties of A.</param>
1079	/// <param name="balance">true: permute A in order to increase the
1080	/// numerical stability, false: do not permute A.</param>
1081	/// <returns>eigenvalues as vector (if V is null) or as diagonoal
1082	/// matrix (if V was requested, i.e. not null).</returns>
1083	/// <remarks><para>The eigenvalues of A are found by use of the
1084	/// Lapack functions dgeevx, sgeevx, cgeevx and zgeevx. </para>
1085	/// <para>The arrays returned will be of complex inner type,
1086	/// since no further constraints are set on the structure of
1087	/// A (it may be nonsymmetric). Use
1088	/// <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>)"/>
1089	/// or <see cref="ILNumerics.ILMath.eigSymm(ILInArray<double>,ILOutArray<double>)"/>
1090	/// functions for computing the real eigenvalues of symmetric
1091	/// matrices explicitly.</para>
1092	/// <para>Depending on the parameter <paramref name="balance"/>,
1093	/// A will be balanced first. This includes permutations and
1094	/// scaling of A in order to improve the conditioning of the
1095	/// eigenvalues.</para></remarks>
1096	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>)"/>
1097	/// <seealso cref="ILNumerics.ILMath.eig(ILInArray<double>,ILOutArray<complex>,ref MatrixProperties,bool)"/>
1098	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if a
1099	/// is not square</exception>
1100	public static ILRetArray< complex > eig(ILInArray< complex > A, ILOutArray< complex > V, ref MatrixProperties propsA, bool balance) {
1101	using (ILScope.Enter(A)) {
1102	if (A.IsEmpty) {
1103	if (!object.Equals(V,null))
1104	V.a = empty<complex>(A.Size);
1105	return empty<complex>(A.Size);
1106	}
1107	ILArray< complex > ret = empty< complex >(ILSize.Empty00);
1108	int n = A.Size[0];
1109	bool createVR = (object.Equals(V,null))? false:true;
1110	if (n != A.Size[1])
1111	throw new ILArgumentException("eig: matrix A must be square");
1112	propsA \|= MatrixProperties.Square;
1113	if (((propsA & MatrixProperties.Hermitian) != 0 \|\| ILMath.ishermitian(A.C))) {
1114	propsA \|= MatrixProperties.Hermitian;
1115	ILArray< complex > Vd = null;
1116	if (createVR)
1117	Vd = returnType< complex >();
1118	ILArray< complex > tmpRet = eigSymm(A,Vd);
1119	if (createVR)
1120	V.a = (Vd);
1121	ret = (tmpRet);
1122	} else {
1123	// nonsymmetric case
1124	char bal = (balance)? 'B':'N', jobvr;
1125	ILRetArray< complex > tmpA = A.C;
1126	complex [] vr = null;
1127	complex [] wr = ILMemoryPool.Pool.New< complex >(n);
1128
1129	double [] scale = ILMemoryPool.Pool.New< double >(n);
1130	double [] rconde = ILMemoryPool.Pool.New< double >(n);
1131	double [] rcondv = ILMemoryPool.Pool.New< double >(n);
1132	double abnrm = 0;
1133	int ldvr, ilo = 0, ihi = 0, info = 0;
1134	if (createVR) {
1135	ldvr = n;
1136	vr = ILMemoryPool.Pool.New< complex >(n * n);
1137	jobvr = 'V';
1138	} else {
1139	ldvr = 1;
1140	vr = new complex [1];
1141	jobvr = 'N';
1142	}
1143	Lapack.zgeevx(bal,'N',jobvr,'N',n,tmpA.GetArrayForWrite(),n,wr, new complex [1],1,vr,ldvr,ref ilo,ref ihi,scale,ref abnrm,rconde,rcondv,ref info);
1144	ILMemoryPool.Pool.Free(rconde);
1145	ILMemoryPool.Pool.Free(rcondv);
1146	ILMemoryPool.Pool.Free(scale);
1147	if (info != 0)
1148	throw new ILArgumentException("eig: error in Lapack '?geevx': (" + info + ")");
1149	// create eigenvalues
1150	complex [] retArr = ILMemoryPool.Pool.New< complex >(n);
1151	for (int i = 0; i < n; i++) {
1152	retArr[i] = wr[i];
1153	}
1154	ret = new ILArray< complex > (retArr,n,1);
1155	if (createVR) {
1156	V.a = array<complex> (vr,n,n);
1157	if (createVR)
1158	ret.a = ILMath.diag< complex>(ret);
1159	}
1160
1161	ILMemoryPool.Pool.Free(vr);
1162	ILMemoryPool.Pool.Free(wr);
1163	}
1164	return ret;
1165	}
1166	}
1167
1168	/// <summary>
1169	/// Find all eigenvalues of symmetric (hermitian) matrix
1170	/// </summary>
1171	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
1172	/// <returns>Array of size [n,1] with eigenvalues of A.</returns>
1173	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
1174	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
1175	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
1176	public static ILRetArray< complex > eigSymm (ILInArray< complex > A) {
1177	using (ILScope.Enter(A)) {
1178	if (A.IsEmpty) {
1179	return empty< complex>(A.Size);
1180	}
1181	int n = A.Size[0];
1182	if (n != A.Size[1])
1183	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
1184	int m = 0,info = 0;
1185	ILArray< double > w = new ILArray< double > (new double [n],1,n);
1186	complex [] z = new complex [1]; ;
1187	int [] isuppz = new int[2 * n];
1188	complex [] AcArr = A.C.GetArrayForWrite();
1189	Lapack.zheevr ('N','A','U',n,AcArr,n,0,0,0,0,0,ref m,w.GetArrayForWrite(),z,1,isuppz,ref info);
1190	ILMemoryPool.Pool.Free(AcArr);
1191	return ILMath.tocomplex(w);
1192	}
1193	}
1194	/// <summary>
1195	/// Find all eigenvalues and -vectors of symmetric (hermitian) matrix
1196	/// </summary>
1197	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
1198	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors
1199	/// will not be computed and V is not changed. In order to make the function return the vectors, V should be initiialized with ILMath.returnType before calling eigSymm.</param>
1200	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
1201	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
1202	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
1203	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square.</exception>
1204	public static ILRetArray< complex > eigSymm (ILInArray< complex > A, ILOutArray< complex > V) {
1205	using (ILScope.Enter(A)) {
1206	if (A.IsEmpty) {
1207	if (!object.Equals(V,null))
1208	V.a = empty<complex>(A.Size);
1209	return empty<complex>(A.Size);
1210	}
1211	int n = A.Size[0];
1212	if (n != A.Size[1])
1213	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
1214	int m = 0,ldz = 0,info = 0;
1215	ILArray< double > w = new ILArray< double >(ILMemoryPool.Pool.New< double >(n),n,1);
1216	complex [] z;
1217	char jobz;
1218	if (object.Equals(V,null)) {
1219	z = new complex [1];
1220	jobz = 'N';
1221	ldz = 1;
1222	} else {
1223	z = ILMemoryPool.Pool.New< complex >(n * n);
1224	jobz = 'V';
1225	ldz = n;
1226	}
1227	int [] isuppz = ILMemoryPool.Pool.New<int>( 2 * n);
1228	complex [] AcArr = A.C.GetArrayForWrite();
1229	Lapack.zheevr (jobz,'A','U',n,AcArr,n,1,n,0,0,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
1230	ILMemoryPool.Pool.Free(AcArr);
1231	ILMemoryPool.Pool.Free(isuppz);
1232	if (info != 0)
1233	throw new ILException("error returned from lapack: " + info);
1234	if (jobz == 'V') {
1235	System.Diagnostics.Debug.Assert(!object.Equals(V,null));
1236	V.a = array< complex > (z,n,n);
1237	V.a = V[full,r(0,m-1)];
1238	return ILMath.diag( ILMath.tocomplex(w));
1239	} else {
1240	ILMemoryPool.Pool.Free(z);
1241	return ILMath.tocomplex(w);
1242	}
1243	}
1244	}
1245	/// <summary>
1246	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
1247	/// </summary>
1248	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
1249	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
1250	/// <param name="rangeStart">Specify the lowest limit for the range of eigenvalues to be queried.</param>
1251	/// <param name="rangeEnd">Specify the upper limit for the range of eigenvalues to be queried.</param>
1252	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
1253	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
1254	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
1255	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/></exception>
1256	public static ILRetArray< complex > eigSymm (ILInArray< complex > A, ILOutArray< complex > V, int rangeStart, int rangeEnd) {
1257	using (ILScope.Enter(A)) {
1258	if (A.IsEmpty) {
1259	if (!object.Equals(V,null))
1260	V.a = empty<complex>(A.Size);
1261	return empty<complex>(A.Size);
1262	}
1263	int n = A.Size[0];
1264	if (n != A.Size[1])
1265	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian.");
1266	int m = 0,ldz = 0,info = 0;
1267	if (rangeEnd < rangeStart \|\| rangeStart < 1)
1268	throw new ILArgumentException("invalid range of eigenvalues requested");
1269	ILArray< double > w = array< double > (
1270	ILMemoryPool.Pool.New< double>(n),1,n);
1271	complex [] z;
1272	char jobz;
1273	if (object.Equals(V,null)) {
1274	z = new complex [1];
1275	jobz = 'N';
1276	ldz = 1;
1277	} else {
1278	z = ILMemoryPool.Pool.New<complex>(n * n);
1279	jobz = 'V';
1280	ldz = n;
1281	}
1282	int [] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
1283	complex [] AcArr = A.C.GetArrayForWrite();
1284	Lapack.zheevr (jobz,'I','U',n,AcArr,n,0,0,rangeStart,rangeEnd,0,ref m,w.GetArrayForWrite(),z,ldz,isuppz,ref info);
1285	ILMemoryPool.Pool.Free(isuppz);
1286	ILMemoryPool.Pool.Free(AcArr);
1287
1288	if (jobz == 'V') {
1289	V.a = array< complex >(z,n,n);
1290	V.a = V[full,r(0,m-1)];
1291	}
1292	ILMemoryPool.Pool.Free(z);
1293	return ILMath.diag( ILMath.tocomplex(w));
1294	}
1295	}
1296	/// <summary>
1297	/// Find some eigenvalues and -vectors of symmetric (hermitian) matrix
1298	/// </summary>
1299	/// <param name="A">Input matrix, Size [n x n], symmetric (hermitian for complex A) </param>
1300	/// <param name="V">Output: n eigenvectors as columns. Size [n x n]. If V is null on input, the eigenvectors will not be computed and V is not changed. </param>
1301	/// <param name="rangeStart">The eigenvalues will be returned by increasing size. This will determine the number of the first eigenvalue to be returned.</param>
1302	/// <param name="rangeEnd">Determine the number of the last eigenvalue to be returned.</param>
1303	/// <returns>Diagonal matrix of size [n,n] with eigenvalues of A on the main diagonal.</returns>
1304	/// <remarks><para>For computation the Lapack functions dsyevr, ssyevr, chesvr and zheesv are used. </para>
1305	/// <para>Since A is symmetric, the eigenvalues will always be real. Therefore the return value will be of the same inner type as A.</para></remarks>
1306	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A is not square or <paramref name="rangeEnd"/> < <paramref name="rangeStart"/> or if either one is <= 0.</exception>
1307	public static ILRetArray< complex> eigSymm( ILInArray< complex> A, ILOutArray< complex> V, double rangeStart, double rangeEnd ) {
1308	using (ILScope.Enter(A)) {
1309	if (A.IsEmpty) {
1310	if (!object.Equals(V, null))
1311	V.a = empty<complex>(A.Size);
1312	return empty<complex>(A.Size);
1313	}
1314	int n = A.Size[0];
1315	if (n != A.Size[1])
1316	throw new ILArgumentException("input matrix A must be square and symmetric/hermitian");
1317	int m = 0, ldz = 0, info = 0;
1318	if (rangeStart > rangeEnd)
1319	throw new ILArgumentException("invalid range of eigenvalues requested");
1320	ILArray< double> w = zeros<double>(1, n);
1321
1322	complex[] z;
1323	char jobz;
1324	if (object.Equals(V, null)) {
1325	z = new complex[1];
1326	jobz = 'N';
1327	ldz = 1;
1328	} else {
1329	z = ILMemoryPool.Pool.New<complex>(n * n);
1330	jobz = 'V';
1331	ldz = n;
1332	}
1333	int[] isuppz = ILMemoryPool.Pool.New<int>(2 * n);
1334
1335
1336	complex[] AcArr = A.C.GetArrayForWrite();
1337
1338	Lapack.zheevr(jobz, 'V', 'U', n, AcArr, n, rangeStart, rangeEnd, 0, 0, 0, ref m, w.GetArrayForWrite(), z, ldz, isuppz, ref info);
1339	ILMemoryPool.Pool.Free(AcArr);
1340	ILMemoryPool.Pool.Free(isuppz);
1341
1342	if (jobz == 'V') {
1343	V.a = array< complex>(z, n, n);
1344	V.a = V[full, r(0, m - 1)];
1345	}
1346	ILMemoryPool.Pool.Free(z);
1347	return diag( ILMath.tocomplex(w));
1348	}
1349	}
1350
1351	#endregion HYCALPER AUTO GENERATED CODE
1352
1353	/// <summary>
1354	/// Specifies the type of eigenproblem
1355	/// </summary>
1356	/// <remarks>The enumeration describes possible problem definitions for generelized eigenproblems:
1357	/// <list type="bullet">
1358	/// <item>Ax_eq_lambBx: AV = rB*V</item>
1359	/// <item>ABx_eq_lambx: ABV = r*V</item>
1360	/// <item>BAx_eq_lambx: BAV = r*V</item>
1361	/// </list></remarks>
1362	public enum GenEigenType {
1363	/// <summary>
1364	/// AV = rB*V
1365	/// </summary>
1366	Ax_eq_lambBx,
1367	/// <summary>
1368	/// ABV = r*V
1369	/// </summary>
1370	ABx_eq_lambx,
1371	/// <summary>
1372	/// BAV = r*V
1373	/// </summary>
1374	BAx_eq_lambx
1375	}
1376
1377
1378
1379
1380	/// <summary>
1381	/// Compute eigenvalues <it>lambda</it> of symmetrical/hermitian inputs A and B: AV=lamdaB*V
1382	/// </summary>
1383	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1384	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1385	/// <returns>Vector of eigenvalues. size [n x 1]</returns>
1386	/// <remarks>
1387	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1388	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1389	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1390	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1391	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if either A and/or B was found not to be symmetric/hermitian</exception>
1392	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1393	public static ILRetArray<double> eigSymm(ILInArray<double> A, ILInArray<double> B) {
1394	using (ILScope.Enter(A, B)) {
1395	return eigSymm(A, B, null, GenEigenType.Ax_eq_lambBx, false);
1396	}
1397	}
1398
1399	/// <summary>
1400	/// Compute eigenvalues and eigenvectors of symmetric/hermitian input
1401	/// </summary>
1402	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1403	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1404	/// <param name="outV">[Output] Returns eigenvectors in columns (size [n x n]). </param>
1405	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1406	/// <returns>Vector of eigenvalues. The return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1407	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1408	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1409	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1410	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1411	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1412	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1413	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1414	public static ILRetArray<double> eigSymm(ILInArray<double> A, ILInArray<double> B, ILOutArray<double> outV, bool skipSymmCheck) {
1415	using (ILScope.Enter(A, B)) {
1416	return eigSymm(A, B, outV, GenEigenType.Ax_eq_lambBx, skipSymmCheck);
1417	}
1418	}
1419
1420	/// <summary>
1421	/// Compute eigenvalues and eigenvectors (optional) of symmetric/hermitian input
1422	/// </summary>
1423	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1424	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1425	/// <param name="outV">[Output] If on input not null-> returns eigenvectors in columns (size [n x n]). If null on input -> eigenvectors will not get computed.</param>
1426	/// <param name="type">Determine the type of problem. This is one of the following types:
1427	/// <list type="bullet">
1428	/// <item>Ax_eq_lambBx: AV = rB*V</item>
1429	/// <item>ABx_eq_lambx: ABV = r*V</item>
1430	/// <item>BAx_eq_lambx: BAV = r*V</item>
1431	/// </list>Here 'r' is the eigenvalue corresponding to the eigenvector 'V'.</param>
1432	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1433	/// <returns>Vector of eigenvalues. If the eigenvectors are requested as well (V not null on input),
1434	/// the return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1435	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1436	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1437	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1438	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1439	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1440	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1441	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1442	public static ILRetArray<double> eigSymm(ILInArray<double> A, ILInArray<double> B, ILOutArray< double> outV, GenEigenType type, bool skipSymmCheck) {
1443	using (ILScope.Enter(A, B)) {
1444	// check input arguments
1445	if (object.Equals(A, null) \|\| object.Equals(B, null))
1446	throw new ILArgumentException("A and B must not be null!");
1447	if (!A.IsMatrix \|\| !B.IsMatrix)
1448	throw new ILArgumentException("A & B must be matrices!");
1449	int n = A.Size[0];
1450	if (n != A.Size[1])
1451	throw new ILArgumentException("input matrices must be square!");
1452	if (!A.Size.IsSameSize(B.Size))
1453	throw new ILArgumentException("A and B must have the same size!");
1454	if (A.IsEmpty) {
1455	if (object.Equals(outV, null))
1456	return empty< double>(A.Size);
1457	else {
1458	outV.a = A.C;
1459	return empty< double>(A.Size);
1460	}
1461	}
1462	if (!skipSymmCheck && !ILMath.ishermitian(A)) {
1463	throw new ILArgumentException("A must be hermitian!");
1464	}
1465	if (!skipSymmCheck && !ILMath.ishermitian(B)) {
1466	throw new ILArgumentException("B must be hermitian!");
1467	}
1468	int info = -1;
1469	int itype = 1;
1470	switch (type) {
1471	case GenEigenType.Ax_eq_lambBx:
1472	itype = 1;
1473	break;
1474	case GenEigenType.ABx_eq_lambx:
1475	itype = 2;
1476	break;
1477	case GenEigenType.BAx_eq_lambx:
1478	itype = 3;
1479	break;
1480	}
1481	char jobz = 'N';
1482	if (!object.Equals(outV, null)) {
1483	jobz = 'V';
1484	outV.a = copyUpperTriangle(A, n, n);
1485	} else {
1486	ILArray< double> tmpOutV = copyUpperTriangle(A, n, n);
1487	outV = tmpOutV;
1488	}
1489	ILArray< double> BC = B.C;
1490
1491	double[] w = ILMemoryPool.Pool.New< double>(n);
1492	/!HC:lapack_func/
1493	Lapack.dsygv(itype, jobz, 'U', n, outV.GetArrayForWrite(), n, BC.GetArrayForWrite(), BC.Size[0], w, ref info);
1494	if (info == 0) {
1495	if (jobz == 'N')
1496	return array< double>(w, 1, n);
1497	else
1498	return diag(array< double>(w, 1, n));
1499	} else if (info < 0) {
1500	throw new ILArgumentException("invalid parameter reported from Lapack module: #" + (-info));
1501	} else {
1502	if (info <= n) {
1503	throw new ILArgumentException(String.Format("eigSymm did not converge! {0} off-diagonal elements unequal 0", info));
1504	} else if (info < 2 * n) {
1505	throw new ILArgumentException("eigSymm: B must be positive definite!");
1506	} else {
1507	throw new ILArgumentException("eigSymm: unknown Lapack module error");
1508	}
1509	}
1510	}
1511	}
1512
1513	#region HYCALPER AUTO GENERATED CODE
1514
1515
1516
1517	/// <summary>
1518	/// Compute eigenvalues <it>lambda</it> of symmetrical/hermitian inputs A and B: AV=lamdaB*V
1519	/// </summary>
1520	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1521	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1522	/// <returns>Vector of eigenvalues. size [n x 1]</returns>
1523	/// <remarks>
1524	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1525	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1526	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1527	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1528	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if either A and/or B was found not to be symmetric/hermitian</exception>
1529	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1530	public static ILRetArray<float> eigSymm(ILInArray<float> A, ILInArray<float> B) {
1531	using (ILScope.Enter(A, B)) {
1532	return eigSymm(A, B, null, GenEigenType.Ax_eq_lambBx, false);
1533	}
1534	}
1535
1536	/// <summary>
1537	/// Compute eigenvalues and eigenvectors of symmetric/hermitian input
1538	/// </summary>
1539	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1540	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1541	/// <param name="outV">[Output] Returns eigenvectors in columns (size [n x n]). </param>
1542	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1543	/// <returns>Vector of eigenvalues. The return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1544	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1545	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1546	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1547	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1548	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1549	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1550	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1551	public static ILRetArray<float> eigSymm(ILInArray<float> A, ILInArray<float> B, ILOutArray<float> outV, bool skipSymmCheck) {
1552	using (ILScope.Enter(A, B)) {
1553	return eigSymm(A, B, outV, GenEigenType.Ax_eq_lambBx, skipSymmCheck);
1554	}
1555	}
1556
1557	/// <summary>
1558	/// Compute eigenvalues and eigenvectors (optional) of symmetric/hermitian input
1559	/// </summary>
1560	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1561	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1562	/// <param name="outV">[Output] If on input not null-> returns eigenvectors in columns (size [n x n]). If null on input -> eigenvectors will not get computed.</param>
1563	/// <param name="type">Determine the type of problem. This is one of the following types:
1564	/// <list type="bullet">
1565	/// <item>Ax_eq_lambBx: AV = rB*V</item>
1566	/// <item>ABx_eq_lambx: ABV = r*V</item>
1567	/// <item>BAx_eq_lambx: BAV = r*V</item>
1568	/// </list>Here 'r' is the eigenvalue corresponding to the eigenvector 'V'.</param>
1569	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1570	/// <returns>Vector of eigenvalues. If the eigenvectors are requested as well (V not null on input),
1571	/// the return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1572	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1573	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1574	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1575	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1576	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1577	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1578	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1579	public static ILRetArray<float> eigSymm(ILInArray<float> A, ILInArray<float> B, ILOutArray< float> outV, GenEigenType type, bool skipSymmCheck) {
1580	using (ILScope.Enter(A, B)) {
1581	// check input arguments
1582	if (object.Equals(A, null) \|\| object.Equals(B, null))
1583	throw new ILArgumentException("A and B must not be null!");
1584	if (!A.IsMatrix \|\| !B.IsMatrix)
1585	throw new ILArgumentException("A & B must be matrices!");
1586	int n = A.Size[0];
1587	if (n != A.Size[1])
1588	throw new ILArgumentException("input matrices must be square!");
1589	if (!A.Size.IsSameSize(B.Size))
1590	throw new ILArgumentException("A and B must have the same size!");
1591	if (A.IsEmpty) {
1592	if (object.Equals(outV, null))
1593	return empty< float>(A.Size);
1594	else {
1595	outV.a = A.C;
1596	return empty< float>(A.Size);
1597	}
1598	}
1599	if (!skipSymmCheck && !ILMath.ishermitian(A)) {
1600	throw new ILArgumentException("A must be hermitian!");
1601	}
1602	if (!skipSymmCheck && !ILMath.ishermitian(B)) {
1603	throw new ILArgumentException("B must be hermitian!");
1604	}
1605	int info = -1;
1606	int itype = 1;
1607	switch (type) {
1608	case GenEigenType.Ax_eq_lambBx:
1609	itype = 1;
1610	break;
1611	case GenEigenType.ABx_eq_lambx:
1612	itype = 2;
1613	break;
1614	case GenEigenType.BAx_eq_lambx:
1615	itype = 3;
1616	break;
1617	}
1618	char jobz = 'N';
1619	if (!object.Equals(outV, null)) {
1620	jobz = 'V';
1621	outV.a = copyUpperTriangle(A, n, n);
1622	} else {
1623	ILArray< float> tmpOutV = copyUpperTriangle(A, n, n);
1624	outV = tmpOutV;
1625	}
1626	ILArray< float> BC = B.C;
1627
1628	float[] w = ILMemoryPool.Pool.New< float>(n);
1629
1630	Lapack.ssygv(itype, jobz, 'U', n, outV.GetArrayForWrite(), n, BC.GetArrayForWrite(), BC.Size[0], w, ref info);
1631	if (info == 0) {
1632	if (jobz == 'N')
1633	return array< float>(w, 1, n);
1634	else
1635	return diag(array< float>(w, 1, n));
1636	} else if (info < 0) {
1637	throw new ILArgumentException("invalid parameter reported from Lapack module: #" + (-info));
1638	} else {
1639	if (info <= n) {
1640	throw new ILArgumentException(String.Format("eigSymm did not converge! {0} off-diagonal elements unequal 0", info));
1641	} else if (info < 2 * n) {
1642	throw new ILArgumentException("eigSymm: B must be positive definite!");
1643	} else {
1644	throw new ILArgumentException("eigSymm: unknown Lapack module error");
1645	}
1646	}
1647	}
1648	}
1649
1650
1651	/// <summary>
1652	/// Compute eigenvalues <it>lambda</it> of symmetrical/hermitian inputs A and B: AV=lamdaB*V
1653	/// </summary>
1654	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1655	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1656	/// <returns>Vector of eigenvalues. size [n x 1]</returns>
1657	/// <remarks>
1658	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1659	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1660	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1661	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1662	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if either A and/or B was found not to be symmetric/hermitian</exception>
1663	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1664	public static ILRetArray<float> eigSymm(ILInArray<fcomplex> A, ILInArray<fcomplex> B) {
1665	using (ILScope.Enter(A, B)) {
1666	return eigSymm(A, B, null, GenEigenType.Ax_eq_lambBx, false);
1667	}
1668	}
1669
1670	/// <summary>
1671	/// Compute eigenvalues and eigenvectors of symmetric/hermitian input
1672	/// </summary>
1673	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1674	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1675	/// <param name="outV">[Output] Returns eigenvectors in columns (size [n x n]). </param>
1676	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1677	/// <returns>Vector of eigenvalues. The return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1678	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1679	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1680	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1681	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1682	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1683	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1684	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1685	public static ILRetArray<float> eigSymm(ILInArray<fcomplex> A, ILInArray<fcomplex> B, ILOutArray<fcomplex> outV, bool skipSymmCheck) {
1686	using (ILScope.Enter(A, B)) {
1687	return eigSymm(A, B, outV, GenEigenType.Ax_eq_lambBx, skipSymmCheck);
1688	}
1689	}
1690
1691	/// <summary>
1692	/// Compute eigenvalues and eigenvectors (optional) of symmetric/hermitian input
1693	/// </summary>
1694	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1695	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1696	/// <param name="outV">[Output] If on input not null-> returns eigenvectors in columns (size [n x n]). If null on input -> eigenvectors will not get computed.</param>
1697	/// <param name="type">Determine the type of problem. This is one of the following types:
1698	/// <list type="bullet">
1699	/// <item>Ax_eq_lambBx: AV = rB*V</item>
1700	/// <item>ABx_eq_lambx: ABV = r*V</item>
1701	/// <item>BAx_eq_lambx: BAV = r*V</item>
1702	/// </list>Here 'r' is the eigenvalue corresponding to the eigenvector 'V'.</param>
1703	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1704	/// <returns>Vector of eigenvalues. If the eigenvectors are requested as well (V not null on input),
1705	/// the return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1706	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1707	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1708	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1709	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1710	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1711	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1712	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1713	public static ILRetArray<float> eigSymm(ILInArray<fcomplex> A, ILInArray<fcomplex> B, ILOutArray< fcomplex> outV, GenEigenType type, bool skipSymmCheck) {
1714	using (ILScope.Enter(A, B)) {
1715	// check input arguments
1716	if (object.Equals(A, null) \|\| object.Equals(B, null))
1717	throw new ILArgumentException("A and B must not be null!");
1718	if (!A.IsMatrix \|\| !B.IsMatrix)
1719	throw new ILArgumentException("A & B must be matrices!");
1720	int n = A.Size[0];
1721	if (n != A.Size[1])
1722	throw new ILArgumentException("input matrices must be square!");
1723	if (!A.Size.IsSameSize(B.Size))
1724	throw new ILArgumentException("A and B must have the same size!");
1725	if (A.IsEmpty) {
1726	if (object.Equals(outV, null))
1727	return empty< float>(A.Size);
1728	else {
1729	outV.a = A.C;
1730	return empty< float>(A.Size);
1731	}
1732	}
1733	if (!skipSymmCheck && !ILMath.ishermitian(A)) {
1734	throw new ILArgumentException("A must be hermitian!");
1735	}
1736	if (!skipSymmCheck && !ILMath.ishermitian(B)) {
1737	throw new ILArgumentException("B must be hermitian!");
1738	}
1739	int info = -1;
1740	int itype = 1;
1741	switch (type) {
1742	case GenEigenType.Ax_eq_lambBx:
1743	itype = 1;
1744	break;
1745	case GenEigenType.ABx_eq_lambx:
1746	itype = 2;
1747	break;
1748	case GenEigenType.BAx_eq_lambx:
1749	itype = 3;
1750	break;
1751	}
1752	char jobz = 'N';
1753	if (!object.Equals(outV, null)) {
1754	jobz = 'V';
1755	outV.a = copyUpperTriangle(A, n, n);
1756	} else {
1757	ILArray< fcomplex> tmpOutV = copyUpperTriangle(A, n, n);
1758	outV = tmpOutV;
1759	}
1760	ILArray< fcomplex> BC = B.C;
1761
1762	float[] w = ILMemoryPool.Pool.New< float>(n);
1763
1764	Lapack.chegv(itype, jobz, 'U', n, outV.GetArrayForWrite(), n, BC.GetArrayForWrite(), BC.Size[0], w, ref info);
1765	if (info == 0) {
1766	if (jobz == 'N')
1767	return array< float>(w, 1, n);
1768	else
1769	return diag(array< float>(w, 1, n));
1770	} else if (info < 0) {
1771	throw new ILArgumentException("invalid parameter reported from Lapack module: #" + (-info));
1772	} else {
1773	if (info <= n) {
1774	throw new ILArgumentException(String.Format("eigSymm did not converge! {0} off-diagonal elements unequal 0", info));
1775	} else if (info < 2 * n) {
1776	throw new ILArgumentException("eigSymm: B must be positive definite!");
1777	} else {
1778	throw new ILArgumentException("eigSymm: unknown Lapack module error");
1779	}
1780	}
1781	}
1782	}
1783
1784
1785	/// <summary>
1786	/// Compute eigenvalues <it>lambda</it> of symmetrical/hermitian inputs A and B: AV=lamdaB*V
1787	/// </summary>
1788	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1789	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1790	/// <returns>Vector of eigenvalues. size [n x 1]</returns>
1791	/// <remarks>
1792	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1793	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1794	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1795	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1796	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if either A and/or B was found not to be symmetric/hermitian</exception>
1797	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1798	public static ILRetArray<double> eigSymm(ILInArray<complex> A, ILInArray<complex> B) {
1799	using (ILScope.Enter(A, B)) {
1800	return eigSymm(A, B, null, GenEigenType.Ax_eq_lambBx, false);
1801	}
1802	}
1803
1804	/// <summary>
1805	/// Compute eigenvalues and eigenvectors of symmetric/hermitian input
1806	/// </summary>
1807	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1808	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1809	/// <param name="outV">[Output] Returns eigenvectors in columns (size [n x n]). </param>
1810	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1811	/// <returns>Vector of eigenvalues. The return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1812	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1813	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1814	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1815	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1816	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1817	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1818	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1819	public static ILRetArray<double> eigSymm(ILInArray<complex> A, ILInArray<complex> B, ILOutArray<complex> outV, bool skipSymmCheck) {
1820	using (ILScope.Enter(A, B)) {
1821	return eigSymm(A, B, outV, GenEigenType.Ax_eq_lambBx, skipSymmCheck);
1822	}
1823	}
1824
1825	/// <summary>
1826	/// Compute eigenvalues and eigenvectors (optional) of symmetric/hermitian input
1827	/// </summary>
1828	/// <param name="A">Square, symmetric/hermitian input matrix, size [n x n]</param>
1829	/// <param name="B">Square, symmetric/hermitian and positive definite matrix, size [n x n]</param>
1830	/// <param name="outV">[Output] If on input not null-> returns eigenvectors in columns (size [n x n]). If null on input -> eigenvectors will not get computed.</param>
1831	/// <param name="type">Determine the type of problem. This is one of the following types:
1832	/// <list type="bullet">
1833	/// <item>Ax_eq_lambBx: AV = rB*V</item>
1834	/// <item>ABx_eq_lambx: ABV = r*V</item>
1835	/// <item>BAx_eq_lambx: BAV = r*V</item>
1836	/// </list>Here 'r' is the eigenvalue corresponding to the eigenvector 'V'.</param>
1837	/// <param name="skipSymmCheck">true: skip tests for A and B being hermitian.</param>
1838	/// <returns>Vector of eigenvalues. If the eigenvectors are requested as well (V not null on input),
1839	/// the return value will be a diagonal matrix with the eigenvalues on the main diagonal.</returns>
1840	/// <remarks><para>The eigenvectors in 'V' are not normalized!</para>
1841	/// <para>Internally, the generalized eigenproblem AV = rBV will be reduced to B<sup>-1</sup>AV = rV using cholesky factorization. The
1842	/// computations are handled by LAPACK functions DSYGV,SSYGV,CHEGV and ZHEGV.</para></remarks>
1843	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if B was not positive definite</exception>
1844	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if A and B was not of the same size</exception>
1845	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if <paramref name="skipSymmCheck"/> is false and either A and/or B was found not to be symmetric/hermitian</exception>
1846	/// <exception cref="ILNumerics.Exceptions.ILArgumentException">if the algorithm did not converge. All exceptions will contain an informational message describing the problem verbosely.</exception>
1847	public static ILRetArray<double> eigSymm(ILInArray<complex> A, ILInArray<complex> B, ILOutArray< complex> outV, GenEigenType type, bool skipSymmCheck) {
1848	using (ILScope.Enter(A, B)) {
1849	// check input arguments
1850	if (object.Equals(A, null) \|\| object.Equals(B, null))
1851	throw new ILArgumentException("A and B must not be null!");
1852	if (!A.IsMatrix \|\| !B.IsMatrix)
1853	throw new ILArgumentException("A & B must be matrices!");
1854	int n = A.Size[0];
1855	if (n != A.Size[1])
1856	throw new ILArgumentException("input matrices must be square!");
1857	if (!A.Size.IsSameSize(B.Size))
1858	throw new ILArgumentException("A and B must have the same size!");
1859	if (A.IsEmpty) {
1860	if (object.Equals(outV, null))
1861	return empty< double>(A.Size);
1862	else {
1863	outV.a = A.C;
1864	return empty< double>(A.Size);
1865	}
1866	}
1867	if (!skipSymmCheck && !ILMath.ishermitian(A)) {
1868	throw new ILArgumentException("A must be hermitian!");
1869	}
1870	if (!skipSymmCheck && !ILMath.ishermitian(B)) {
1871	throw new ILArgumentException("B must be hermitian!");
1872	}
1873	int info = -1;
1874	int itype = 1;
1875	switch (type) {
1876	case GenEigenType.Ax_eq_lambBx:
1877	itype = 1;
1878	break;
1879	case GenEigenType.ABx_eq_lambx:
1880	itype = 2;
1881	break;
1882	case GenEigenType.BAx_eq_lambx:
1883	itype = 3;
1884	break;
1885	}
1886	char jobz = 'N';
1887	if (!object.Equals(outV, null)) {
1888	jobz = 'V';
1889	outV.a = copyUpperTriangle(A, n, n);
1890	} else {
1891	ILArray< complex> tmpOutV = copyUpperTriangle(A, n, n);
1892	outV = tmpOutV;
1893	}
1894	ILArray< complex> BC = B.C;
1895
1896	double[] w = ILMemoryPool.Pool.New< double>(n);
1897
1898	Lapack.zhegv(itype, jobz, 'U', n, outV.GetArrayForWrite(), n, BC.GetArrayForWrite(), BC.Size[0], w, ref info);
1899	if (info == 0) {
1900	if (jobz == 'N')
1901	return array< double>(w, 1, n);
1902	else
1903	return diag(array< double>(w, 1, n));
1904	} else if (info < 0) {
1905	throw new ILArgumentException("invalid parameter reported from Lapack module: #" + (-info));
1906	} else {
1907	if (info <= n) {
1908	throw new ILArgumentException(String.Format("eigSymm did not converge! {0} off-diagonal elements unequal 0", info));
1909	} else if (info < 2 * n) {
1910	throw new ILArgumentException("eigSymm: B must be positive definite!");
1911	} else {
1912	throw new ILArgumentException("eigSymm: unknown Lapack module error");
1913	}
1914	}
1915	}
1916	}
1917
1918	#endregion HYCALPER AUTO GENERATED CODE
1919	}
1920	}
1921

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

Download in other formats:

Update cookies preferences