Lapack_Eigenvalues.cxx
1 // Copyright (C) 2003-2009 Marc Duruflé
2 // Copyright (C) 2001-2009 Vivien Mallet
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_LAPACK_EIGENVALUES_CXX
22 
23 
24 #include "Lapack_Eigenvalues.hxx"
25 
26 
27 namespace Seldon
28 {
29 
30 
32  // STANDARD EIGENVALUE PROBLEM //
33 
34 
35  /* RowMajor */
36 
37 
38  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
39  void GetEigenvalues(Matrix<float, Prop, RowMajor, Allocator1>& A,
40  Vector<float, VectFull, Allocator2>& wr,
41  Vector<float, VectFull, Allocator3>& wi,
42  LapackInfo& info)
43  {
44 #ifdef SELDON_CHECK_DIMENSIONS
45  if (A.GetM() != A.GetN())
46  throw WrongDim("GetEigenvalues", "Matrix must be squared");
47 #endif
48 
49  int n = A.GetM(), lwork = 6*n;
50  char jobvl('N');
51  char jobvr('N');
52  Vector<float> work(lwork);
53  wr.Reallocate(n);
54  wi.Reallocate(n);
55  sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
56  A.GetData(), &n, A.GetData(), &n, work.GetData(),
57  &lwork, &info.GetInfoRef());
58 
59 #ifdef SELDON_LAPACK_CHECK_INFO
60  if (info.GetInfo() != 0)
61  throw LapackError(info.GetInfo(), "GetEigenvalues",
62  "Failed to find eigenvalues ");
63 #endif
64 
65  }
66 
67 
68  template<class Prop, class Allocator1, class Allocator2,
69  class Allocator3, class Allocator4>
70  void GetEigenvaluesEigenvectors(Matrix<float, Prop, RowMajor, Allocator1>& A,
71  Vector<float, VectFull, Allocator2>& wr,
72  Vector<float, VectFull, Allocator3>& wi,
73  Matrix<float, General, RowMajor, Allocator4>& zr,
74  LapackInfo& info)
75  {
76 #ifdef SELDON_CHECK_DIMENSIONS
77  if (A.GetM() != A.GetN())
78  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
79 #endif
80 
81  int n = A.GetM(), lwork = 6*n;
82  char jobvl('V');
83  char jobvr('N');
84  Vector<float> work(lwork);
85  wr.Reallocate(n);
86  wi.Reallocate(n);
87  zr.Reallocate(n, n);
88  sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
89  zr.GetData(), &n, zr.GetData(), &n, work.GetData(), &lwork,
90  &info.GetInfoRef());
91 
92 #ifdef SELDON_LAPACK_CHECK_INFO
93  if (info.GetInfo() != 0)
94  throw LapackError(info.GetInfo(), "GetEigenvalues",
95  "Failed to find eigenvalues ");
96 #endif
97 
98  Transpose(zr);
99  // conjugate if necessary
100  int i = 0;
101  while (i < n)
102  {
103  if (i < (n-1))
104  if (wr(i) == wr(i+1))
105  {
106  for (int j = 0; j < n; j++)
107  zr(j,i+1) = -zr(j,i+1);
108 
109  i++;
110  }
111 
112  i++;
113  }
114  }
115 
116 
117  template<class Prop, class Allocator1, class Allocator2>
118  void GetEigenvalues(Matrix<complex<float>, Prop, RowMajor, Allocator1>& A,
119  Vector<complex<float>, VectFull, Allocator2>& w,
120  LapackInfo& info)
121  {
122 #ifdef SELDON_CHECK_DIMENSIONS
123  if (A.GetM() != A.GetN())
124  throw WrongDim("GetEigenvalues", "Matrix must be squared");
125 #endif
126 
127  int n = A.GetM();
128  char jobl('N'), jobr('N'); int lwork = 3*n;
129  Vector<complex<float> > work(lwork);
130  Vector<float> rwork(2*n);
131  w.Reallocate(n);
132  cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
133  A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), &lwork,
134  rwork.GetData(), &info.GetInfoRef());
135 
136 #ifdef SELDON_LAPACK_CHECK_INFO
137  if (info.GetInfo() != 0)
138  throw LapackError(info.GetInfo(), "GetEigenvalues",
139  "Failed to find eigenvalues ");
140 #endif
141 
142  }
143 
144 
145  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
146  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
147  Prop, RowMajor, Allocator1>& A,
148  Vector<complex<float>, VectFull, Allocator2>& w,
149  Matrix<complex<float>,
150  General, RowMajor, Allocator3>& z,
151  LapackInfo& info)
152  {
153 #ifdef SELDON_CHECK_DIMENSIONS
154  if (A.GetM() != A.GetN())
155  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
156 #endif
157 
158  int n = A.GetM();
159  char jobl('V'), jobr('N'); int lwork = 3*n;
160  Vector<complex<float> > work(lwork);
161  Vector<float> rwork(2*n);
162  w.Reallocate(n);
163  z.Reallocate(n, n);
164  cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
165  z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(), &lwork,
166  rwork.GetData(), &info.GetInfoRef());
167 
168 #ifdef SELDON_LAPACK_CHECK_INFO
169  if (info.GetInfo() != 0)
170  throw LapackError(info.GetInfo(), "GetEigenvalues",
171  "Failed to find eigenvalues ");
172 #endif
173 
174  TransposeConj(z);
175  }
176 
177 
178  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
179  void GetEigenvalues(Matrix<double, Prop, RowMajor, Allocator1>& A,
180  Vector<double, VectFull, Allocator2>& wr,
181  Vector<double, VectFull, Allocator3>& wi,
182  LapackInfo& info)
183  {
184 #ifdef SELDON_CHECK_DIMENSIONS
185  if (A.GetM() != A.GetN())
186  throw WrongDim("GetEigenvalues", "Matrix must be squared");
187 #endif
188 
189  int n = A.GetM(), lwork = 6*n;
190  char jobvl('N'), jobvr('N');
191  Vector<double> work(lwork);
192  wr.Reallocate(n);
193  wi.Reallocate(n);
194  dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
195  A.GetData(), &n, A.GetData(), &n, work.GetData(),
196  &lwork, &info.GetInfoRef());
197 
198 #ifdef SELDON_LAPACK_CHECK_INFO
199  if (info.GetInfo() != 0)
200  throw LapackError(info.GetInfo(), "GetEigenvalues",
201  "Failed to find eigenvalues ");
202 #endif
203 
204  }
205 
206 
207  template<class Prop, class Allocator1, class Allocator2,
208  class Allocator3, class Allocator4>
209  void GetEigenvaluesEigenvectors(Matrix<double, Prop, RowMajor, Allocator1>& A,
210  Vector<double, VectFull, Allocator2>& wr,
211  Vector<double, VectFull, Allocator3>& wi,
212  Matrix<double, General, RowMajor, Allocator4>& zr,
213  LapackInfo& info)
214  {
215 #ifdef SELDON_CHECK_DIMENSIONS
216  if (A.GetM() != A.GetN())
217  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
218 #endif
219 
220  int n = A.GetM(), lwork = 6*n;
221  char jobvl('V'), jobvr('N');
222  Vector<double> work(lwork);
223  wr.Reallocate(n);
224  wi.Reallocate(n);
225  zr.Reallocate(n, n);
226  dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
227  zr.GetData(), &n, zr.GetData(), &n, work.GetData(), &lwork,
228  &info.GetInfoRef());
229 
230 #ifdef SELDON_LAPACK_CHECK_INFO
231  if (info.GetInfo() != 0)
232  throw LapackError(info.GetInfo(), "GetEigenvalues",
233  "Failed to find eigenvalues ");
234 #endif
235 
236  Transpose(zr);
237  // conjugate if necessary
238  int i = 0;
239  while (i < n)
240  {
241  if (i < (n-1))
242  if (wr(i) == wr(i+1))
243  {
244  for (int j = 0; j < n; j++)
245  zr(j,i+1) = -zr(j,i+1);
246 
247  i++;
248  }
249 
250  i++;
251  }
252  }
253 
254 
255  template<class Prop, class Allocator1, class Allocator2>
256  void GetEigenvalues(Matrix<complex<double>, Prop, RowMajor, Allocator1>& A,
257  Vector<complex<double>, VectFull, Allocator2>& w,
258  LapackInfo& info)
259  {
260 #ifdef SELDON_CHECK_DIMENSIONS
261  if (A.GetM() != A.GetN())
262  throw WrongDim("GetEigenvalues", "Matrix must be squared");
263 #endif
264 
265  int n = A.GetM();
266  char jobl('N'), jobr('N'); int lwork = 3*n;
267  Vector<complex<double> > work(lwork);
268  Vector<double> rwork(2*n);
269  w.Reallocate(n);
270  zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
271  A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(), &lwork,
272  rwork.GetData(), &info.GetInfoRef());
273 
274 #ifdef SELDON_LAPACK_CHECK_INFO
275  if (info.GetInfo() != 0)
276  throw LapackError(info.GetInfo(), "GetEigenvalues",
277  "Failed to find eigenvalues ");
278 #endif
279 
280  }
281 
282 
283  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
284  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
285  Prop, RowMajor, Allocator1>& A,
286  Vector<complex<double>,
287  VectFull, Allocator2>& w,
288  Matrix<complex<double>,
289  General, RowMajor, Allocator3>& z,
290  LapackInfo& info)
291  {
292 #ifdef SELDON_CHECK_DIMENSIONS
293  if (A.GetM() != A.GetN())
294  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
295 #endif
296 
297  int n = A.GetM();
298  char jobl('V'), jobr('N'); int lwork = 3*n;
299  Vector<complex<double> > work(lwork);
300  Vector<double> rwork(2*n);
301  w.Reallocate(n);
302  z.Reallocate(n, n);
303  zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
304  z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(),
305  &lwork, rwork.GetData(), &info.GetInfoRef());
306 
307 #ifdef SELDON_LAPACK_CHECK_INFO
308  if (info.GetInfo() != 0)
309  throw LapackError(info.GetInfo(), "GetEigenvalues",
310  "Failed to find eigenvalues ");
311 #endif
312 
313  TransposeConj(z);
314  }
315 
316 
317  /* ColMajor */
318 
319 
320  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
321  void GetEigenvalues(Matrix<float, Prop, ColMajor, Allocator1>& A,
322  Vector<float, VectFull, Allocator2>& wr,
323  Vector<float, VectFull, Allocator3>& wi,
324  LapackInfo& info)
325  {
326 #ifdef SELDON_CHECK_DIMENSIONS
327  if (A.GetM() != A.GetN())
328  throw WrongDim("GetEigenvalues", "Matrix must be squared");
329 #endif
330 
331  int n = A.GetM(), lwork = 6*n;
332  char jobvl('N'), jobvr('N');
333  Vector<float> work(lwork);
334  wr.Reallocate(n);
335  wi.Reallocate(n);
336  sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
337  A.GetData(), &n, A.GetData(), &n, work.GetData(), &lwork,
338  &info.GetInfoRef());
339 
340 #ifdef SELDON_LAPACK_CHECK_INFO
341  if (info.GetInfo() != 0)
342  throw LapackError(info.GetInfo(), "GetEigenvalues",
343  "Failed to find eigenvalues ");
344 #endif
345 
346  }
347 
348  template<class Prop, class Allocator1, class Allocator2,
349  class Allocator3, class Allocator4>
350  void GetEigenvaluesEigenvectors(Matrix<float, Prop, ColMajor, Allocator1>& A,
351  Vector<float, VectFull, Allocator2>& wr,
352  Vector<float, VectFull, Allocator3>& wi,
353  Matrix<float, General, ColMajor, Allocator4>&zr,
354  LapackInfo& info)
355  {
356 #ifdef SELDON_CHECK_DIMENSIONS
357  if (A.GetM() != A.GetN())
358  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
359 #endif
360 
361  int n = A.GetM(), lwork = 6*n;
362  char jobvl('N'), jobvr('V');
363  Vector<float> work(lwork);
364  wr.Reallocate(n);
365  wi.Reallocate(n);
366  zr.Reallocate(n, n);
367  sgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
368  zr.GetData(), &n, zr.GetData(), &n, work.GetData(),
369  &lwork, &info.GetInfoRef());
370 
371 #ifdef SELDON_LAPACK_CHECK_INFO
372  if (info.GetInfo() != 0)
373  throw LapackError(info.GetInfo(), "GetEigenvalues",
374  "Failed to find eigenvalues ");
375 #endif
376 
377  }
378 
379 
380  template<class Prop, class Allocator1, class Allocator2>
381  void GetEigenvalues(Matrix<complex<float>, Prop, ColMajor, Allocator1>& A,
382  Vector<complex<float>, VectFull, Allocator2>& w,
383  LapackInfo& info)
384  {
385 #ifdef SELDON_CHECK_DIMENSIONS
386  if (A.GetM() != A.GetN())
387  throw WrongDim("GetEigenvalues", "Matrix must be squared");
388 #endif
389 
390  int n = A.GetM();
391  char jobl('N'), jobr('N'); int lwork = 3*n;
392  Vector<complex<float> > work(lwork);
393  Vector<float> rwork(2*n);
394  w.Reallocate(n);
395  cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
396  A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(),
397  &lwork, rwork.GetData(), &info.GetInfoRef());
398 
399 #ifdef SELDON_LAPACK_CHECK_INFO
400  if (info.GetInfo() != 0)
401  throw LapackError(info.GetInfo(), "GetEigenvalues",
402  "Failed to find eigenvalues ");
403 #endif
404 
405  }
406 
407 
408  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
409  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
410  Prop, ColMajor, Allocator1>& A,
411  Vector<complex<float>,
412  VectFull, Allocator2>& w,
413  Matrix<complex<float>,
414  General, ColMajor, Allocator3>& z,
415  LapackInfo& info)
416  {
417 #ifdef SELDON_CHECK_DIMENSIONS
418  if (A.GetM() != A.GetN())
419  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
420 #endif
421 
422  int n = A.GetM();
423  char jobl('N'), jobr('V'); int lwork = 3*n;
424  Vector<complex<float> > work(lwork);
425  Vector<float> rwork(2*n);
426  w.Reallocate(n);
427  z.Reallocate(n, n);
428  cgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
429  z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(),
430  &lwork, rwork.GetData(), &info.GetInfoRef());
431 
432 #ifdef SELDON_LAPACK_CHECK_INFO
433  if (info.GetInfo() != 0)
434  throw LapackError(info.GetInfo(), "GetEigenvalues",
435  "Failed to find eigenvalues ");
436 #endif
437 
438  }
439 
440  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
441  void GetEigenvalues(Matrix<double, Prop, ColMajor, Allocator1>& A,
442  Vector<double, VectFull, Allocator2>& wr,
443  Vector<double, VectFull, Allocator3>& wi,
444  LapackInfo& info)
445  {
446 #ifdef SELDON_CHECK_DIMENSIONS
447  if (A.GetM() != A.GetN())
448  throw WrongDim("GetEigenvalues", "Matrix must be squared");
449 #endif
450 
451  int n = A.GetM(), lwork = 6*n;
452  char jobvl('N'), jobvr('N');
453  Vector<double> work(lwork);
454  wr.Reallocate(n);
455  wi.Reallocate(n);
456  dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
457  A.GetData(), &n, A.GetData(), &n, work.GetData(), &lwork,
458  &info.GetInfoRef());
459 
460 #ifdef SELDON_LAPACK_CHECK_INFO
461  if (info.GetInfo() != 0)
462  throw LapackError(info.GetInfo(), "GetEigenvalues",
463  "Failed to find eigenvalues ");
464 #endif
465 
466  }
467 
468 
469  template<class Prop, class Allocator1, class Allocator2,
470  class Allocator3, class Allocator4>
471  void GetEigenvaluesEigenvectors(Matrix<double, Prop, ColMajor, Allocator1>& A,
472  Vector<double, VectFull, Allocator2>& wr,
473  Vector<double, VectFull, Allocator3>& wi,
474  Matrix<double, General, ColMajor, Allocator4>&zr,
475  LapackInfo& info)
476  {
477 #ifdef SELDON_CHECK_DIMENSIONS
478  if (A.GetM() != A.GetN())
479  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
480 #endif
481 
482  int n = A.GetM(), lwork = 6*n;
483  char jobvl('N'), jobvr('V');
484  Vector<double> work(lwork);
485  wr.Reallocate(n);
486  wi.Reallocate(n);
487  zr.Reallocate(n, n);
488  dgeev_(&jobvl, &jobvr, &n, A.GetData(), &n, wr.GetData(), wi.GetData(),
489  zr.GetData(), &n, zr.GetData(), &n, work.GetData(),
490  &lwork, &info.GetInfoRef());
491 
492 #ifdef SELDON_LAPACK_CHECK_INFO
493  if (info.GetInfo() != 0)
494  throw LapackError(info.GetInfo(), "GetEigenvalues",
495  "Failed to find eigenvalues ");
496 #endif
497 
498  }
499 
500 
501  template<class Prop, class Allocator1, class Allocator2>
502  void GetEigenvalues(Matrix<complex<double>, Prop, ColMajor, Allocator1>& A,
503  Vector<complex<double>, VectFull, Allocator2>& w,
504  LapackInfo& info)
505  {
506 #ifdef SELDON_CHECK_DIMENSIONS
507  if (A.GetM() != A.GetN())
508  throw WrongDim("GetEigenvalues", "Matrix must be squared");
509 #endif
510 
511  int n = A.GetM();
512  char jobl('N'), jobr('N'); int lwork = 3*n;
513  Vector<complex<double> > work(lwork);
514  Vector<double> rwork(2*n);
515  w.Reallocate(n);
516  zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
517  A.GetDataVoid(), &n, A.GetData(), &n, work.GetDataVoid(),
518  &lwork, rwork.GetData(), &info.GetInfoRef());
519 
520 #ifdef SELDON_LAPACK_CHECK_INFO
521  if (info.GetInfo() != 0)
522  throw LapackError(info.GetInfo(), "GetEigenvalues",
523  "Failed to find eigenvalues ");
524 #endif
525 
526  }
527 
528 
529  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
530  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
531  Prop, ColMajor, Allocator1>& A,
532  Vector<complex<double>,
533  VectFull, Allocator2>& w,
534  Matrix<complex<double>,
535  General, ColMajor, Allocator3>& z,
536  LapackInfo& info)
537  {
538 #ifdef SELDON_CHECK_DIMENSIONS
539  if (A.GetM() != A.GetN())
540  throw WrongDim("GetEigenvaluesEigenvectors", "Matrix must be squared");
541 #endif
542 
543  int n = A.GetM();
544  char jobl('N'), jobr('V'); int lwork = 3*n;
545  Vector<complex<double> > work(lwork);
546  Vector<double> rwork(2*n);
547  w.Reallocate(n);
548  z.Reallocate(n, n);
549  zgeev_(&jobl, &jobr, &n, A.GetDataVoid(), &n, w.GetDataVoid(),
550  z.GetDataVoid(), &n, z.GetData(), &n, work.GetDataVoid(),
551  &lwork, rwork.GetData(), &info.GetInfoRef());
552 
553 #ifdef SELDON_LAPACK_CHECK_INFO
554  if (info.GetInfo() != 0)
555  throw LapackError(info.GetInfo(), "GetEigenvalues",
556  "Failed to find eigenvalues ");
557 #endif
558 
559  }
560 
561 
562  /* RowSym */
563 
564 
565  template<class Prop, class Allocator1, class Allocator2>
566  void GetEigenvalues(Matrix<float, Prop, RowSym, Allocator1>& A,
567  Vector<float, VectFull, Allocator2>& w,
568  LapackInfo& info)
569  {
570  int n = A.GetM();
571  char uplo('L'); char job('N');
572  int lwork = 3*n;
573  Vector<float> work(lwork);
574  w.Reallocate(n);
575  ssyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(),
576  &lwork, &info.GetInfoRef());
577 
578 #ifdef SELDON_LAPACK_CHECK_INFO
579  if (info.GetInfo() != 0)
580  throw LapackError(info.GetInfo(), "GetEigenvalues",
581  "Failed to find eigenvalues ");
582 #endif
583 
584  }
585 
586 
587  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
588  void GetEigenvaluesEigenvectors(Matrix<float, Prop, RowSym, Allocator1>& A,
589  Vector<float, VectFull, Allocator2>& w,
590  Matrix<float, General, RowMajor, Allocator3>& z,
591  LapackInfo& info)
592  {
593  int n = A.GetM();
594  char uplo('L'); char job('V');
595  int lwork = 3*n;
596  z.Reallocate(n, n);
597  for (int i = 0; i < n; i++)
598  for (int j = 0; j < n; j++)
599  z(i,j) = A(i,j);
600 
601  A.Clear();
602  Vector<float> work(lwork);
603  w.Reallocate(n);
604  ssyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), work.GetData(),
605  &lwork, &info.GetInfoRef());
606 
607 #ifdef SELDON_LAPACK_CHECK_INFO
608  if (info.GetInfo() != 0)
609  throw LapackError(info.GetInfo(), "GetEigenvalues",
610  "Failed to find eigenvalues ");
611 #endif
612 
613  Transpose(z);
614  }
615 
616  template<class Prop, class Allocator1, class Allocator2>
617  void GetEigenvalues(Matrix<complex<float>, Prop, RowSym, Allocator1>& A,
618  Vector<complex<float>, VectFull, Allocator2>& w,
619  LapackInfo& info)
620  {
621  int n = A.GetM();
622  Matrix<complex<float>, General, ColMajor> B(n,n);
623  for (int i = 0; i < n; i++)
624  for (int j = 0; j < n; j++)
625  B(i,j) = A(i,j);
626 
627  A.Clear();
628  GetEigenvalues(B, w, info);
629  }
630 
631  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
632  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
633  Prop, RowSym, Allocator1>& A,
634  Vector<complex<float>,
635  VectFull, Allocator2>& w,
636  Matrix<complex<float>,
637  General, RowMajor, Allocator3>& z,
638  LapackInfo& info)
639  {
640  int n = A.GetM();
641  Matrix<complex<float>, General, RowMajor> B(n,n);
642  for (int i = 0; i < n; i++)
643  for (int j = 0; j < n; j++)
644  B(i,j) = A(i,j);
645 
646  A.Clear();
647  GetEigenvaluesEigenvectors(B, w, z, info);
648  }
649 
650 
651  template<class Prop, class Allocator1, class Allocator2>
652  void GetEigenvalues(Matrix<double, Prop, RowSym, Allocator1>& A,
653  Vector<double, VectFull, Allocator2>& w,
654  LapackInfo& info)
655  {
656  int n = A.GetM();
657  char uplo('L'), job('N');
658  int lwork = 3*n;
659  Vector<double> work(lwork);
660  w.Reallocate(n);
661  dsyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(),
662  &lwork, &info.GetInfoRef());
663 
664 #ifdef SELDON_LAPACK_CHECK_INFO
665  if (info.GetInfo() != 0)
666  throw LapackError(info.GetInfo(), "GetEigenvalues",
667  "Failed to find eigenvalues ");
668 #endif
669 
670  }
671 
672 
673  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
674  void GetEigenvaluesEigenvectors(Matrix<double, Prop, RowSym, Allocator1>& A,
675  Vector<double, VectFull, Allocator2>& w,
676  Matrix<double, General, RowMajor, Allocator3>& z,
677  LapackInfo& info)
678  {
679  int n = A.GetM();
680  char uplo('L'); char job('V');
681  int lwork = 3*n;
682  z.Reallocate(n, n);
683  for (int i = 0; i < n; i++)
684  for (int j = 0; j < n; j++)
685  z(i,j) = A(i,j);
686 
687  A.Clear();
688  Vector<double> work(lwork);
689  w.Reallocate(n);
690  dsyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(), work.GetData(),
691  &lwork, &info.GetInfoRef());
692 
693 #ifdef SELDON_LAPACK_CHECK_INFO
694  if (info.GetInfo() != 0)
695  throw LapackError(info.GetInfo(), "GetEigenvalues",
696  "Failed to find eigenvalues ");
697 #endif
698 
699  Transpose(z);
700  }
701 
702  template<class Prop, class Allocator1, class Allocator2>
703  void GetEigenvalues(Matrix<complex<double>, Prop, RowSym, Allocator1>& A,
704  Vector<complex<double>, VectFull, Allocator2>& w,
705  LapackInfo& info)
706  {
707  int n = A.GetM();
708  Matrix<complex<double>, General, ColMajor> B(n,n);
709  for (int i = 0; i < n; i++)
710  for (int j = 0; j < n; j++)
711  B(i,j) = A(i,j);
712 
713  A.Clear();
714  GetEigenvalues(B, w, info);
715  }
716 
717  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
718  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
719  Prop, RowSym, Allocator1>& A,
720  Vector<complex<double>,
721  VectFull, Allocator2>& w,
722  Matrix<complex<double>,
723  General, RowMajor, Allocator3>& z,
724  LapackInfo& info)
725  {
726  int n = A.GetM();
727  Matrix<complex<double>, General, RowMajor> B(n,n);
728  for (int i = 0; i < n; i++)
729  for (int j = 0; j < n; j++)
730  B(i,j) = A(i,j);
731 
732  A.Clear();
733  w.Reallocate(n);
734  z.Reallocate(n, n);
735  GetEigenvaluesEigenvectors(B, w, z, info);
736  }
737 
738 
739  /* ColSym */
740 
741 
742  template<class Prop, class Allocator1, class Allocator2>
743  void GetEigenvalues(Matrix<float, Prop, ColSym, Allocator1>& A,
744  Vector<float, VectFull, Allocator2>& w,
745  LapackInfo& info)
746  {
747  int n = A.GetM();
748  char uplo('U'); char job('N');
749  int lwork = 3*n;
750  Vector<float> work(lwork);
751  w.Reallocate(n);
752  ssyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(),
753  &lwork, &info.GetInfoRef());
754 
755 #ifdef SELDON_LAPACK_CHECK_INFO
756  if (info.GetInfo() != 0)
757  throw LapackError(info.GetInfo(), "GetEigenvalues",
758  "Failed to find eigenvalues ");
759 #endif
760 
761  }
762 
763 
764  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
765  void GetEigenvaluesEigenvectors(Matrix<float, Prop, ColSym, Allocator1>& A,
766  Vector<float, VectFull, Allocator2>& w,
767  Matrix<float, General, ColMajor, Allocator3>& z,
768  LapackInfo& info)
769  {
770  int n = A.GetM();
771  char uplo('U'); char job('V');
772  int lwork = 3*n;
773  z.Reallocate(n, n);
774  for (int i = 0; i < n; i++)
775  for (int j = 0; j < n; j++)
776  z(i,j) = A(i,j);
777 
778  A.Clear();
779  Vector<float> work(lwork);
780  w.Reallocate(n);
781  ssyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(),
782  work.GetData(), &lwork, &info.GetInfoRef());
783 
784 #ifdef SELDON_LAPACK_CHECK_INFO
785  if (info.GetInfo() != 0)
786  throw LapackError(info.GetInfo(), "GetEigenvalues",
787  "Failed to find eigenvalues ");
788 #endif
789 
790  }
791 
792 
793  template<class Prop, class Allocator1, class Allocator2>
794  void GetEigenvalues(Matrix<complex<float>, Prop, ColSym, Allocator1>& A,
795  Vector<complex<float>, VectFull, Allocator2>& w,
796  LapackInfo& info)
797  {
798  int n = A.GetM();
799  Matrix<complex<float>, General, ColMajor> B(n,n);
800  for (int i = 0; i < n; i++)
801  for (int j = 0; j < n; j++)
802  B(i,j) = A(i,j);
803 
804  A.Clear();
805  GetEigenvalues(B, w, info);
806  }
807 
808 
809  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
810  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
811  Prop, ColSym, Allocator1>& A,
812  Vector<complex<float>,
813  VectFull, Allocator2>& w,
814  Matrix<complex<float>,
815  General, ColMajor, Allocator3>& z,
816  LapackInfo& info)
817  {
818  int n = A.GetM();
819  Matrix<complex<float>, General, ColMajor> B(n,n);
820  for (int i = 0; i < n; i++)
821  for (int j = 0; j < n; j++)
822  B(i,j) = A(i,j);
823 
824  A.Clear();
825  w.Reallocate(n);
826  z.Reallocate(n, n);
827  GetEigenvaluesEigenvectors(B, w, z, info);
828  }
829 
830 
831  template<class Prop, class Allocator1, class Allocator2>
832  void GetEigenvalues(Matrix<double, Prop, ColSym, Allocator1>& A,
833  Vector<double, VectFull, Allocator2>& w,
834  LapackInfo& info)
835  {
836  int n = A.GetM();
837  char uplo('U'), job('N');
838  int lwork = 3*n;
839  Vector<double> work(lwork);
840  w.Reallocate(n);
841  dsyev_(&job, &uplo, &n, A.GetData(), &n, w.GetData(), work.GetData(),
842  &lwork, &info.GetInfoRef());
843 
844 #ifdef SELDON_LAPACK_CHECK_INFO
845  if (info.GetInfo() != 0)
846  throw LapackError(info.GetInfo(), "GetEigenvalues",
847  "Failed to find eigenvalues ");
848 #endif
849 
850  }
851 
852 
853  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
854  void GetEigenvaluesEigenvectors(Matrix<double, Prop, ColSym, Allocator1>& A,
855  Vector<double, VectFull, Allocator2>& w,
856  Matrix<double, General, ColMajor, Allocator3>& z,
857  LapackInfo& info)
858  {
859  int n = A.GetM();
860  char uplo('U'), job('V');
861  int lwork = 3*n;
862  z.Reallocate(n, n);
863  for (int i = 0; i < n; i++)
864  for (int j = 0; j < n; j++)
865  z(i,j) = A(i,j);
866 
867  A.Clear();
868  Vector<double> work(lwork);
869  w.Reallocate(n);
870 
871  dsyev_(&job, &uplo, &n, z.GetData(), &n, w.GetData(),
872  work.GetData(), &lwork, &info.GetInfoRef());
873 
874 #ifdef SELDON_LAPACK_CHECK_INFO
875  if (info.GetInfo() != 0)
876  throw LapackError(info.GetInfo(), "GetEigenvalues",
877  "Failed to find eigenvalues ");
878 #endif
879 
880  }
881 
882 
883  template<class Prop, class Allocator1, class Allocator2>
884  void GetEigenvalues(Matrix<complex<double>, Prop, ColSym, Allocator1>& A,
885  Vector<complex<double>, VectFull, Allocator2>& w,
886  LapackInfo& info)
887  {
888  int n = A.GetM();
889  Matrix<complex<double>, General, ColMajor> B(n,n);
890  for (int i = 0; i < n; i++)
891  for (int j = 0; j < n; j++)
892  B(i,j) = A(i,j);
893 
894  A.Clear();
895  w.Reallocate(n);
896  GetEigenvalues(B, w, info);
897  }
898 
899 
900  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
901  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
902  Prop, ColSym, Allocator1>& A,
903  Vector<complex<double>,
904  VectFull, Allocator2>& w,
905  Matrix<complex<double>,
906  General, ColMajor, Allocator3>& z,
907  LapackInfo& info)
908  {
909  int n = A.GetM();
910  Matrix<complex<double>, General, ColMajor> B(n,n);
911  for (int i = 0; i < n; i++)
912  for (int j = 0; j < n; j++)
913  B(i,j) = A(i,j);
914 
915  A.Clear();
916  w.Reallocate(n);
917  z.Reallocate(n, n);
918  GetEigenvaluesEigenvectors(B, w, z, info);
919  }
920 
921 
922  /* RowHerm */
923 
924 
925  template<class Prop, class Allocator1, class Allocator2>
926  void GetEigenvalues(Matrix<complex<float>, Prop, RowHerm, Allocator1>& A,
927  Vector<float, VectFull, Allocator2>& w,
928  LapackInfo& info)
929  {
930  int n = A.GetM();
931  char uplo('L'), job('N');
932  int lwork = 2*n;
933  Vector<complex<float> > work(lwork);
934  Vector<float> rwork(3*n);
935  w.Reallocate(n);
936  Conjugate(A);
937  cheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(),
938  work.GetDataVoid(), &lwork, rwork.GetData(),
939  &info.GetInfoRef());
940 
941 #ifdef SELDON_LAPACK_CHECK_INFO
942  if (info.GetInfo() != 0)
943  throw LapackError(info.GetInfo(), "GetEigenvalues",
944  "Failed to find eigenvalues ");
945 #endif
946 
947  }
948 
949 
950  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
951  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
952  Prop, RowHerm, Allocator1>& A,
953  Vector<float, VectFull, Allocator2>& w,
954  Matrix<complex<float>,
955  General, RowMajor, Allocator3>& z,
956  LapackInfo& info)
957  {
958  int n = A.GetM();
959  char uplo('L'), job('V');
960  int lwork = 2*n;
961  z.Reallocate(n, n);
962  for (int i = 0; i < n; i++)
963  for (int j = 0; j < n; j++)
964  z(i, j) = conj(A(i, j));
965 
966  A.Clear();
967  w.Reallocate(n);
968  Vector<complex<float> > work(lwork);
969  Vector<float> rwork(3*n);
970 
971  cheev_(&job, &uplo, &n, z.GetDataVoid(), &n,
972  w.GetData(), work.GetDataVoid(),
973  &lwork, rwork.GetData(), &info.GetInfoRef());
974 
975 #ifdef SELDON_LAPACK_CHECK_INFO
976  if (info.GetInfo() != 0)
977  throw LapackError(info.GetInfo(), "GetEigenvalues",
978  "Failed to find eigenvalues ");
979 #endif
980 
981  Transpose(z);
982  }
983 
984  template<class Prop, class Allocator1, class Allocator2>
985  void GetEigenvalues(Matrix<complex<double>, Prop, RowHerm, Allocator1>& A,
986  Vector<double, VectFull, Allocator2>& w,
987  LapackInfo& info)
988  {
989  int n = A.GetM();
990  char uplo('L'), job('N');
991  int lwork = 2*n;
992  Vector<complex<double> > work(lwork);
993  Vector<double> rwork(3*n);
994  w.Reallocate(n);
995  Conjugate(A);
996  zheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(),
997  work.GetDataVoid(), &lwork, rwork.GetData(),
998  &info.GetInfoRef());
999 
1000 #ifdef SELDON_LAPACK_CHECK_INFO
1001  if (info.GetInfo() != 0)
1002  throw LapackError(info.GetInfo(), "GetEigenvalues",
1003  "Failed to find eigenvalues ");
1004 #endif
1005 
1006  }
1007 
1008 
1009  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1010  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1011  Prop, RowHerm, Allocator1>& A,
1012  Vector<double, VectFull, Allocator2>& w,
1013  Matrix<complex<double>,
1014  General, RowMajor, Allocator3>& z,
1015  LapackInfo& info)
1016  {
1017  int n = A.GetM();
1018  char uplo('L'), job('V');
1019  int lwork = 2*n;
1020  z.Reallocate(n, n);
1021  for (int i = 0; i < n; i++)
1022  for (int j = 0; j < n; j++)
1023  z(i,j) = conj(A(i,j));
1024 
1025  A.Clear();
1026  Vector<complex<double> > work(lwork);
1027  Vector<double> rwork(3*n);
1028  w.Reallocate(n);
1029 
1030  zheev_(&job, &uplo,&n, z.GetDataVoid(),&n, w.GetData(), work.GetDataVoid(),
1031  &lwork, rwork.GetData() , &info.GetInfoRef());
1032 
1033 #ifdef SELDON_LAPACK_CHECK_INFO
1034  if (info.GetInfo() != 0)
1035  throw LapackError(info.GetInfo(), "GetEigenvalues",
1036  "Failed to find eigenvalues ");
1037 #endif
1038 
1039  Transpose(z);
1040  }
1041 
1042 
1043  /* ColHerm */
1044 
1045 
1046  template<class Prop, class Allocator1, class Allocator2>
1047  void GetEigenvalues(Matrix<complex<float>, Prop, ColHerm, Allocator1>& A,
1048  Vector<float, VectFull, Allocator2>& w,
1049  LapackInfo& info)
1050  {
1051  int n = A.GetM();
1052  char uplo('U'); char job('N'); int lwork = 2*n;
1053  Vector<complex<float> > work(lwork);
1054  Vector<float> rwork(3*n);
1055  w.Reallocate(n);
1056  cheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(),
1057  work.GetDataVoid(), &lwork, rwork.GetData(),
1058  &info.GetInfoRef());
1059 
1060 #ifdef SELDON_LAPACK_CHECK_INFO
1061  if (info.GetInfo() != 0)
1062  throw LapackError(info.GetInfo(), "GetEigenvalues",
1063  "Failed to find eigenvalues ");
1064 #endif
1065 
1066  }
1067 
1068 
1069  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1070  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
1071  Prop, ColHerm, Allocator1>& A,
1072  Vector<float, VectFull, Allocator2>& w,
1073  Matrix<complex<float>,
1074  General, ColMajor, Allocator3>& z,
1075  LapackInfo& info)
1076  {
1077  int n = A.GetM();
1078  char uplo('U'), job('V');
1079  int lwork = 2*n;
1080  z.Reallocate(n,n);
1081  for (int i = 0; i < n; i++)
1082  for (int j = 0; j < n; j++)
1083  z(i,j) = A(i,j);
1084 
1085  A.Clear();
1086  Vector<complex<float> > work(lwork);
1087  Vector<float> rwork(3*n);
1088  w.Reallocate(n);
1089 
1090  cheev_(&job, &uplo, &n, z.GetDataVoid(), &n,
1091  w.GetData(), work.GetDataVoid(),
1092  &lwork, rwork.GetData() , &info.GetInfoRef());
1093 
1094 #ifdef SELDON_LAPACK_CHECK_INFO
1095  if (info.GetInfo() != 0)
1096  throw LapackError(info.GetInfo(), "GetEigenvalues",
1097  "Failed to find eigenvalues ");
1098 #endif
1099 
1100  }
1101 
1102 
1103  template<class Prop, class Allocator1, class Allocator2>
1104  void GetEigenvalues(Matrix<complex<double>, Prop, ColHerm, Allocator1>& A,
1105  Vector<double, VectFull, Allocator2>& w,
1106  LapackInfo& info)
1107  {
1108  int n = A.GetM();
1109  char uplo('U'), job('N');
1110  int lwork = 2*n;
1111  Vector<complex<double> > work(lwork);
1112  Vector<double> rwork(3*n);
1113  w.Reallocate(n);
1114  zheev_(&job, &uplo, &n, A.GetDataVoid(), &n, w.GetData(),
1115  work.GetDataVoid(), &lwork, rwork.GetData(),
1116  &info.GetInfoRef());
1117 
1118 #ifdef SELDON_LAPACK_CHECK_INFO
1119  if (info.GetInfo() != 0)
1120  throw LapackError(info.GetInfo(), "GetEigenvalues",
1121  "Failed to find eigenvalues ");
1122 #endif
1123 
1124  }
1125 
1126 
1127  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1128  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1129  Prop, ColHerm, Allocator1>& A,
1130  Vector<double, VectFull, Allocator2>& w,
1131  Matrix<complex<double>,
1132  General, ColMajor, Allocator3>& z,
1133  LapackInfo& info)
1134  {
1135  int n = A.GetM();
1136  char uplo('U'), job('V');
1137  int lwork = 2*n;
1138  z.Reallocate(n,n);
1139  for (int i = 0; i < n; i++)
1140  for (int j = 0; j < n; j++)
1141  z(i,j) = A(i,j);
1142 
1143  A.Clear();
1144  Vector<complex<double> > work(lwork);
1145  Vector<double> rwork(3*n);
1146  w.Reallocate(n);
1147 
1148  zheev_(&job, &uplo, &n, z.GetDataVoid(), &n,
1149  w.GetData(), work.GetDataVoid(),
1150  &lwork, rwork.GetData() , &info.GetInfoRef());
1151 
1152 #ifdef SELDON_LAPACK_CHECK_INFO
1153  if (info.GetInfo() != 0)
1154  throw LapackError(info.GetInfo(), "GetEigenvalues",
1155  "Failed to find eigenvalues ");
1156 #endif
1157 
1158  }
1159 
1160 
1161  /* RowSymPacked */
1162 
1163 
1164  template<class Prop, class Allocator1, class Allocator2>
1165  void GetEigenvalues(Matrix<float, Prop, RowSymPacked, Allocator1>& A,
1166  Vector<float, VectFull, Allocator2>& w,
1167  LapackInfo& info)
1168  {
1169  int n = A.GetM();
1170  char uplo('L'); char job('N');
1171  Vector<float> work(3*n);
1172  w.Reallocate(n);
1173  sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), &n,
1174  work.GetData() , &info.GetInfoRef());
1175 
1176 #ifdef SELDON_LAPACK_CHECK_INFO
1177  if (info.GetInfo() != 0)
1178  throw LapackError(info.GetInfo(), "GetEigenvalues",
1179  "Failed to find eigenvalues ");
1180 #endif
1181 
1182  }
1183 
1184 
1185  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1186  void GetEigenvaluesEigenvectors(Matrix<float, Prop,RowSymPacked, Allocator1>& A,
1187  Vector<float, VectFull, Allocator2>& w,
1188  Matrix<float, General, RowMajor, Allocator3>&z,
1189  LapackInfo& info)
1190  {
1191  int n = A.GetM();
1192  char uplo('L'); char job('V');
1193  Vector<float> work(3*n);
1194  w.Reallocate(n);
1195  z.Reallocate(n,n);
1196  sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), &n,
1197  work.GetData() , &info.GetInfoRef());
1198 
1199 #ifdef SELDON_LAPACK_CHECK_INFO
1200  if (info.GetInfo() != 0)
1201  throw LapackError(info.GetInfo(), "GetEigenvalues",
1202  "Failed to find eigenvalues ");
1203 #endif
1204 
1205  Transpose(z);
1206  }
1207 
1208 
1209  template<class Prop, class Allocator1, class Allocator2>
1210  void GetEigenvalues(Matrix<complex<float>, Prop, RowSymPacked, Allocator1>& A,
1211  Vector<complex<float>, VectFull, Allocator2>& w,
1212  LapackInfo& info)
1213  {
1214  int n = A.GetM();
1215  Matrix<complex<float>, General, ColMajor> B(n,n);
1216  for (int i = 0; i < n; i++)
1217  for (int j = 0; j < n; j++)
1218  B(i,j) = A(i,j);
1219 
1220  A.Clear();
1221  w.Reallocate(n);
1222  GetEigenvalues(B, w, info);
1223  }
1224 
1225  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1226  void GetEigenvaluesEigenvectors(Matrix<complex<float>, Prop,
1227  RowSymPacked, Allocator1>& A,
1228  Vector<complex<float>,VectFull, Allocator2>& w,
1229  Matrix<complex<float>,
1230  General, RowMajor, Allocator3>& z,
1231  LapackInfo& info)
1232  {
1233  int n = A.GetM();
1234  Matrix<complex<float>, General, RowMajor> B(n,n);
1235  for (int i = 0; i < n; i++)
1236  for (int j = 0; j < n; j++)
1237  B(i,j) = A(i,j);
1238 
1239  A.Clear();
1240  w.Reallocate(n);
1241  z.Reallocate(n,n);
1242  GetEigenvaluesEigenvectors(B, w, z, info);
1243  }
1244 
1245 
1246  template<class Prop, class Allocator1, class Allocator2>
1247  void GetEigenvalues(Matrix<double, Prop, RowSymPacked, Allocator1>& A,
1248  Vector<double, VectFull, Allocator2>& w,
1249  LapackInfo& info)
1250  {
1251  int n = A.GetM();
1252  char uplo('L'), job('N');
1253  Vector<double> work(3*n);
1254  w.Reallocate(n);
1255  dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(), &n,
1256  work.GetData() , &info.GetInfoRef());
1257 
1258 #ifdef SELDON_LAPACK_CHECK_INFO
1259  if (info.GetInfo() != 0)
1260  throw LapackError(info.GetInfo(), "GetEigenvalues",
1261  "Failed to find eigenvalues ");
1262 #endif
1263 
1264  }
1265 
1266 
1267  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1268  void GetEigenvaluesEigenvectors(Matrix<double,Prop,RowSymPacked, Allocator1>& A,
1269  Vector<double, VectFull, Allocator2>& w,
1270  Matrix<double, General, RowMajor, Allocator3>&z,
1271  LapackInfo& info)
1272  {
1273  int n = A.GetM();
1274  char uplo('L'), job('V');
1275  Vector<double> work(3*n);
1276  w.Reallocate(n);
1277  z.Reallocate(n,n);
1278  dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(), &n,
1279  work.GetData() , &info.GetInfoRef());
1280 
1281 #ifdef SELDON_LAPACK_CHECK_INFO
1282  if (info.GetInfo() != 0)
1283  throw LapackError(info.GetInfo(), "GetEigenvalues",
1284  "Failed to find eigenvalues ");
1285 #endif
1286 
1287  Transpose(z);
1288  }
1289 
1290 
1291  template<class Prop, class Allocator1, class Allocator2>
1292  void GetEigenvalues(Matrix<complex<double>, Prop, RowSymPacked, Allocator1>& A,
1293  Vector<complex<double>, VectFull, Allocator2>& w,
1294  LapackInfo& info)
1295  {
1296  int n = A.GetM();
1297  Matrix<complex<double>, General, ColMajor> B(n,n);
1298  for (int i = 0; i < n; i++)
1299  for (int j = 0; j < n; j++)
1300  B(i,j) = A(i,j);
1301 
1302  A.Clear();
1303  w.Reallocate(n);
1304  GetEigenvalues(B, w, info);
1305  }
1306 
1307 
1308  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1309  void GetEigenvaluesEigenvectors(Matrix<complex<double>, Prop,
1310  RowSymPacked, Allocator1>& A,
1311  Vector<complex<double>,VectFull, Allocator2>& w,
1312  Matrix<complex<double>,
1313  General, RowMajor, Allocator3>& z,
1314  LapackInfo& info)
1315  {
1316  int n = A.GetM();
1317  Matrix<complex<double>, General, RowMajor> B(n,n);
1318  for (int i = 0; i < n; i++)
1319  for (int j = 0; j < n; j++)
1320  B(i,j) = A(i,j);
1321 
1322  A.Clear();
1323  w.Reallocate(n);
1324  z.Reallocate(n,n);
1325  GetEigenvaluesEigenvectors(B, w, z, info);
1326  }
1327 
1328 
1329  /* ColSymPacked */
1330 
1331 
1332  template<class Prop, class Allocator1, class Allocator2>
1333  void GetEigenvalues(Matrix<float, Prop, ColSymPacked, Allocator1>& A,
1334  Vector<float, VectFull, Allocator2>& w,
1335  LapackInfo& info)
1336  {
1337  int n = A.GetM();
1338  char uplo('U'), job('N');
1339  Vector<float> work(3*n);
1340  w.Reallocate(n);
1341  sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(),
1342  &n, work.GetData() , &info.GetInfoRef());
1343 
1344 #ifdef SELDON_LAPACK_CHECK_INFO
1345  if (info.GetInfo() != 0)
1346  throw LapackError(info.GetInfo(), "GetEigenvalues",
1347  "Failed to find eigenvalues ");
1348 #endif
1349 
1350  }
1351 
1352 
1353  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1354  void GetEigenvaluesEigenvectors(Matrix<float, Prop,ColSymPacked, Allocator1>& A,
1355  Vector<float, VectFull, Allocator2>& w,
1356  Matrix<float, General, ColMajor, Allocator3>&z,
1357  LapackInfo& info)
1358  {
1359  int n = A.GetM();
1360  char uplo('U'), job('V');
1361  Vector<float> work(3*n);
1362  w.Reallocate(n);
1363  z.Reallocate(n,n);
1364  sspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(),
1365  &n, work.GetData() , &info.GetInfoRef());
1366 
1367 #ifdef SELDON_LAPACK_CHECK_INFO
1368  if (info.GetInfo() != 0)
1369  throw LapackError(info.GetInfo(), "GetEigenvalues",
1370  "Failed to find eigenvalues ");
1371 #endif
1372 
1373  }
1374 
1375 
1376  template<class Prop, class Allocator1, class Allocator2>
1377  void GetEigenvalues(Matrix<complex<float>,
1378  Prop, ColSymPacked, Allocator1>& A,
1379  Vector<complex<float>, VectFull, Allocator2>& w,
1380  LapackInfo& info)
1381  {
1382  int n = A.GetM();
1383  Matrix<complex<float>, General, ColMajor> B(n,n);
1384  for (int i = 0; i < n; i++)
1385  for (int j = 0; j < n; j++)
1386  B(i,j) = A(i,j);
1387 
1388  A.Clear();
1389  w.Reallocate(n);
1390  GetEigenvalues(B, w, info);
1391  }
1392 
1393 
1394  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1395  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
1396  Prop, ColSymPacked, Allocator1>& A,
1397  Vector<complex<float>, VectFull, Allocator2>& w,
1398  Matrix<complex<float>,
1399  General, ColMajor, Allocator3>& z,
1400  LapackInfo& info)
1401  {
1402  int n = A.GetM();
1403  Matrix<complex<float>, General, ColMajor> B(n,n);
1404  for (int i = 0; i < n; i++)
1405  for (int j = 0; j < n; j++)
1406  B(i,j) = A(i,j);
1407 
1408  A.Clear();
1409  w.Reallocate(n);
1410  z.Reallocate(n,n);
1411  GetEigenvaluesEigenvectors(B, w, z, info);
1412  }
1413 
1414 
1415  template<class Prop, class Allocator1, class Allocator2>
1416  void GetEigenvalues(Matrix<double, Prop, ColSymPacked, Allocator1>& A,
1417  Vector<double, VectFull, Allocator2>& w,
1418  LapackInfo& info)
1419  {
1420  int n = A.GetM();
1421  char uplo('U'), job('N');
1422  Vector<double> work(3*n);
1423  w.Reallocate(n);
1424  dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), A.GetData(),
1425  &n, work.GetData() , &info.GetInfoRef());
1426 
1427 #ifdef SELDON_LAPACK_CHECK_INFO
1428  if (info.GetInfo() != 0)
1429  throw LapackError(info.GetInfo(), "GetEigenvalues",
1430  "Failed to find eigenvalues ");
1431 #endif
1432 
1433  }
1434 
1435 
1436  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1437  void GetEigenvaluesEigenvectors(Matrix<double, Prop, ColSymPacked, Allocator1>& A,
1438  Vector<double, VectFull, Allocator2>& w,
1439  Matrix<double, General, ColMajor, Allocator3>&z,
1440  LapackInfo& info)
1441  {
1442  int n = A.GetM();
1443  char uplo('U'), job('V');
1444  Vector<double> work(3*n);
1445  w.Reallocate(n);
1446  z.Reallocate(n,n);
1447  dspev_(&job, &uplo, &n, A.GetData(), w.GetData(), z.GetData(),
1448  &n, work.GetData() , &info.GetInfoRef());
1449 
1450 #ifdef SELDON_LAPACK_CHECK_INFO
1451  if (info.GetInfo() != 0)
1452  throw LapackError(info.GetInfo(), "GetEigenvalues",
1453  "Failed to find eigenvalues ");
1454 #endif
1455 
1456  }
1457 
1458 
1459  template<class Prop, class Allocator1, class Allocator2>
1460  void GetEigenvalues(Matrix<complex<double>,
1461  Prop, ColSymPacked, Allocator1>& A,
1462  Vector<complex<double>, VectFull, Allocator2>& w,
1463  LapackInfo& info)
1464  {
1465  int n = A.GetM();
1466  Matrix<complex<double>, General, ColMajor> B(n,n);
1467  for (int i = 0; i < n; i++)
1468  for (int j = 0; j < n; j++)
1469  B(i,j) = A(i,j);
1470 
1471  A.Clear();
1472  w.Reallocate(n);
1473  GetEigenvalues(B, w, info);
1474  }
1475 
1476 
1477  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1478  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1479  Prop, ColSymPacked, Allocator1>& A,
1480  Vector<complex<double>, VectFull, Allocator2>& w,
1481  Matrix<complex<double>,
1482  General, ColMajor, Allocator3>& z,
1483  LapackInfo& info)
1484  {
1485  int n = A.GetM();
1486  Matrix<complex<double>, General, ColMajor> B(n,n);
1487  for (int i = 0; i < n; i++)
1488  for (int j = 0; j < n; j++)
1489  B(i,j) = A(i,j);
1490 
1491  A.Clear();
1492  w.Reallocate(n);
1493  z.Reallocate(n,n);
1494  GetEigenvaluesEigenvectors(B, w, z, info);
1495  }
1496 
1497 
1498  /* RowHermPacked */
1499 
1500 
1501  template<class Prop, class Allocator1, class Allocator2>
1502  void GetEigenvalues(Matrix<complex<float>,
1503  Prop, RowHermPacked, Allocator1>& A,
1504  Vector<float, VectFull, Allocator2>& w,
1505  LapackInfo& info)
1506  {
1507  int n = A.GetM();
1508  char uplo('L'), job('N');
1509  Vector<complex<float> > work(2*n);
1510  Vector<float> rwork(3*n);
1511  w.Reallocate(n);
1512  Conjugate(A);
1513  chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n,
1514  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1515 
1516 #ifdef SELDON_LAPACK_CHECK_INFO
1517  if (info.GetInfo() != 0)
1518  throw LapackError(info.GetInfo(), "GetEigenvalues",
1519  "Failed to find eigenvalues ");
1520 #endif
1521 
1522  }
1523 
1524 
1525  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1526  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
1527  Prop, RowHermPacked, Allocator1>& A,
1528  Vector<float, VectFull, Allocator2>& w,
1529  Matrix<complex<float>,
1530  General, RowMajor, Allocator3>&z,
1531  LapackInfo& info)
1532  {
1533  int n = A.GetM();
1534  char uplo('L'), job('V');
1535  Vector<complex<float> > work(2*n);
1536  Vector<float> rwork(3*n);
1537  w.Reallocate(n);
1538  z.Reallocate(n,n);
1539  Conjugate(A);
1540  chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(),
1541  &n, work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1542 
1543 #ifdef SELDON_LAPACK_CHECK_INFO
1544  if (info.GetInfo() != 0)
1545  throw LapackError(info.GetInfo(), "GetEigenvalues",
1546  "Failed to find eigenvalues ");
1547 #endif
1548 
1549  Transpose(z);
1550  }
1551 
1552 
1553  template<class Prop, class Allocator1, class Allocator2>
1554  void GetEigenvalues(Matrix<complex<double>,
1555  Prop, RowHermPacked, Allocator1>& A,
1556  Vector<double, VectFull, Allocator2>& w,
1557  LapackInfo& info)
1558  {
1559  int n = A.GetM();
1560  char uplo('L'), job('N');
1561  Vector<complex<double> > work(2*n);
1562  Vector<double> rwork(3*n);
1563  w.Reallocate(n);
1564  Conjugate(A);
1565  zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n,
1566  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1567 
1568 #ifdef SELDON_LAPACK_CHECK_INFO
1569  if (info.GetInfo() != 0)
1570  throw LapackError(info.GetInfo(), "GetEigenvalues",
1571  "Failed to find eigenvalues ");
1572 #endif
1573 
1574  }
1575 
1576 
1577  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1578  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1579  Prop, RowHermPacked, Allocator1>& A,
1580  Vector<double, VectFull, Allocator2>& w,
1581  Matrix<complex<double>,
1582  General, RowMajor, Allocator3>&z,
1583  LapackInfo& info)
1584  {
1585  int n = A.GetM();
1586  char uplo('L'), job('V');
1587  Vector<complex<double> > work(2*n);
1588  Vector<double> rwork(3*n);
1589  w.Reallocate(n);
1590  z.Reallocate(n, n);
1591  Conjugate(A);
1592  zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(),
1593  &n, work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1594 
1595 #ifdef SELDON_LAPACK_CHECK_INFO
1596  if (info.GetInfo() != 0)
1597  throw LapackError(info.GetInfo(), "GetEigenvalues",
1598  "Failed to find eigenvalues ");
1599 #endif
1600 
1601  Transpose(z);
1602  }
1603 
1604 
1605  /* ColHermPacked */
1606 
1607 
1608  template<class Prop, class Allocator1, class Allocator2>
1609  void GetEigenvalues(Matrix<complex<float>,
1610  Prop, ColHermPacked, Allocator1>& A,
1611  Vector<float, VectFull, Allocator2>& w,
1612  LapackInfo& info)
1613  {
1614  int n = A.GetM();
1615  char uplo('U');
1616  char job('N');
1617  Vector<complex<float> > work(2*n);
1618  Vector<float> rwork(3*n);
1619  w.Reallocate(n);
1620  chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n,
1621  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1622 
1623 #ifdef SELDON_LAPACK_CHECK_INFO
1624  if (info.GetInfo() != 0)
1625  throw LapackError(info.GetInfo(), "GetEigenvalues",
1626  "Failed to find eigenvalues ");
1627 #endif
1628 
1629  }
1630 
1631 
1632  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1633  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
1634  Prop, ColHermPacked, Allocator1>& A,
1635  Vector<float, VectFull, Allocator2>& w,
1636  Matrix<complex<float>,
1637  General, ColMajor, Allocator3>&z,
1638  LapackInfo& info)
1639  {
1640  int n = A.GetM();
1641  char uplo('U'); char job('V');
1642  Vector<complex<float> > work(2*n);
1643  Vector<float> rwork(3*n);
1644  w.Reallocate(n);
1645  z.Reallocate(n,n);
1646  chpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), &n,
1647  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1648 
1649 #ifdef SELDON_LAPACK_CHECK_INFO
1650  if (info.GetInfo() != 0)
1651  throw LapackError(info.GetInfo(), "GetEigenvalues",
1652  "Failed to find eigenvalues ");
1653 #endif
1654 
1655  }
1656 
1657 
1658  template<class Prop, class Allocator1, class Allocator2>
1659  void GetEigenvalues(Matrix<complex<double>,
1660  Prop, ColHermPacked, Allocator1>& A,
1661  Vector<double, VectFull, Allocator2>& w,
1662  LapackInfo& info)
1663  {
1664  int n = A.GetM();
1665  char uplo('U');
1666  char job('N');
1667  Vector<complex<double> > work(2*n);
1668  Vector<double> rwork(3*n);
1669  w.Reallocate(n);
1670  zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), A.GetDataVoid(), &n,
1671  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1672 
1673 #ifdef SELDON_LAPACK_CHECK_INFO
1674  if (info.GetInfo() != 0)
1675  throw LapackError(info.GetInfo(), "GetEigenvalues",
1676  "Failed to find eigenvalues ");
1677 #endif
1678 
1679  }
1680 
1681 
1682  template<class Prop, class Allocator1, class Allocator2, class Allocator3>
1683  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1684  Prop, ColHermPacked, Allocator1>& A,
1685  Vector<double, VectFull, Allocator2>& w,
1686  Matrix<complex<double>,
1687  General, ColMajor, Allocator3>&z,
1688  LapackInfo& info)
1689  {
1690  int n = A.GetM();
1691  char uplo('U'); char job('V');
1692  Vector<complex<double> > work(2*n);
1693  Vector<double> rwork(3*n);
1694  w.Reallocate(n);
1695  z.Reallocate(n,n);
1696  zhpev_(&job, &uplo, &n, A.GetDataVoid(), w.GetData(), z.GetDataVoid(), &n,
1697  work.GetDataVoid(), rwork.GetData(), &info.GetInfoRef());
1698 
1699 #ifdef SELDON_LAPACK_CHECK_INFO
1700  if (info.GetInfo() != 0)
1701  throw LapackError(info.GetInfo(), "GetEigenvalues",
1702  "Failed to find eigenvalues ");
1703 #endif
1704 
1705  }
1706 
1707 
1708  // STANDARD EIGENVALUE PROBLEM //
1710 
1711 
1713  // GENERALIZED EIGENVALUE PROBLEM //
1714 
1715 
1716  /* RowSym */
1717 
1718 
1719  template<class Prop1, class Prop2, class Allocator1,
1720  class Allocator2, class Allocator3>
1721  void GetEigenvalues(Matrix<float, Prop1, RowSym, Allocator1>& A,
1722  Matrix<float, Prop2, RowSym, Allocator2>& B,
1723  Vector<float, VectFull, Allocator3>& w,
1724  LapackInfo& info)
1725  {
1726 #ifdef SELDON_CHECK_DIMENSIONS
1727  if (A.GetM() != B.GetM())
1728  throw WrongDim("GetEigenvalues",
1729  "Matrix A and B must have the same size");
1730 #endif
1731 
1732 
1733  int itype = 1;
1734  int n = A.GetM();
1735  char uplo('L'), job('N');
1736  int lwork = 3*n;
1737  Vector<float> work(lwork);
1738  w.Reallocate(n);
1739  ssygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
1740  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
1741 
1742 #ifdef SELDON_LAPACK_CHECK_INFO
1743  if (info.GetInfo() != 0)
1744  throw LapackError(info.GetInfo(), "GetEigenvalues",
1745  "Failed to find eigenvalues ");
1746 #endif
1747 
1748  }
1749 
1750 
1751  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
1752  class Allocator3, class Allocator4>
1753  void GetEigenvaluesEigenvectors(Matrix<float, Prop1, RowSym, Allocator1>& A,
1754  Matrix<float, Prop2, RowSym, Allocator2>& B,
1755  Vector<float, VectFull, Allocator3>& w,
1756  Matrix<float, General, RowMajor, Allocator4>& z,
1757  LapackInfo& info)
1758  {
1759 #ifdef SELDON_CHECK_DIMENSIONS
1760  if (A.GetM() != B.GetM())
1761  throw WrongDim("GetEigenvaluesEigenvectors",
1762  "Matrix A and B must have the same size");
1763 #endif
1764 
1765  int itype = 1;
1766  int n = A.GetM();
1767  char uplo('L'), job('V');
1768  int lwork = 3*n;
1769  z.Reallocate(n, n);
1770  for (int i = 0; i < n; i++)
1771  for (int j = 0; j < n; j++)
1772  z(i,j) = A(i,j);
1773 
1774  A.Clear();
1775  Vector<float> work(lwork);
1776  w.Reallocate(n);
1777 
1778  ssygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
1779  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
1780 
1781 #ifdef SELDON_LAPACK_CHECK_INFO
1782  if (info.GetInfo() != 0)
1783  throw LapackError(info.GetInfo(), "GetEigenvalues",
1784  "Failed to find eigenvalues ");
1785 #endif
1786 
1787  Transpose(z);
1788  }
1789 
1790 
1791  template<class Prop1, class Prop2, class Allocator1,
1792  class Allocator2, class Allocator4, class Allocator5>
1793  void GetEigenvalues(Matrix<complex<float>, Prop1, RowSym, Allocator1>& A,
1794  Matrix<complex<float>, Prop2, RowSym, Allocator2>& B,
1795  Vector<complex<float>, VectFull, Allocator4>& alpha,
1796  Vector<complex<float>, VectFull, Allocator5>& beta,
1797  LapackInfo& info)
1798  {
1799 #ifdef SELDON_CHECK_DIMENSIONS
1800  if (A.GetM() != B.GetM())
1801  throw WrongDim("GetEigenvalues",
1802  "Matrix A and B must have the same size");
1803 #endif
1804 
1805  int n = A.GetM();
1806  Matrix<complex<float>, General, RowMajor, Allocator1> A2(n, n), B2(n, n);
1807  for (int i = 0; i < n; i++)
1808  for (int j = 0; j < n; j++)
1809  {
1810  A2(i, j) = A(i, j);
1811  B2(i, j) = B(i, j);
1812  }
1813 
1814  A.Clear();
1815  B.Clear();
1816  GetEigenvalues(A2, B2, alpha, beta, info);
1817  }
1818 
1819 
1820  template<class Prop1, class Prop2, class Prop3, class Allocator1,
1821  class Allocator2, class Allocator4,
1822  class Allocator5, class Allocator6>
1823  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
1824  Prop1, RowSym, Allocator1>& A,
1825  Matrix<complex<float>,
1826  Prop2, RowSym, Allocator2>& B,
1827  Vector<complex<float>,
1828  VectFull, Allocator4>& alpha,
1829  Vector<complex<float>,
1830  VectFull, Allocator5>& beta,
1831  Matrix<complex<float>,
1832  Prop3, RowMajor, Allocator6>& V,
1833  LapackInfo& info)
1834  {
1835 #ifdef SELDON_CHECK_DIMENSIONS
1836  if (A.GetM() != B.GetM())
1837  throw WrongDim("GetEigenvaluesEigenvectors",
1838  "Matrix A and B must have the same size");
1839 #endif
1840 
1841  int n = A.GetM();
1842  Matrix<complex<float>, General, RowMajor, Allocator1> A2(n, n), B2(n, n);
1843  for (int i = 0; i < n; i++)
1844  for (int j = 0; j < n; j++)
1845  {
1846  A2(i, j) = A(i, j);
1847  B2(i, j) = B(i, j);
1848  }
1849 
1850  A.Clear();
1851  B.Clear();
1852  GetEigenvaluesEigenvectors(A2, B2, alpha, beta, V, info);
1853  }
1854 
1855 
1856  template<class Prop1, class Prop2, class Allocator1,
1857  class Allocator2, class Allocator3>
1858  void GetEigenvalues(Matrix<double, Prop1, RowSym, Allocator1>& A,
1859  Matrix<double, Prop2, RowSym, Allocator2>& B,
1860  Vector<double, VectFull, Allocator3>& w,
1861  LapackInfo& info)
1862  {
1863 #ifdef SELDON_CHECK_DIMENSIONS
1864  if (A.GetM() != B.GetM())
1865  throw WrongDim("GetEigenvalues",
1866  "Matrix A and B must have the same size");
1867 #endif
1868 
1869  int itype = 1;
1870  int n = A.GetM();
1871  char uplo('L'), job('N');
1872  int lwork = 3*n;
1873  Vector<double> work(lwork);
1874  w.Reallocate(n);
1875  dsygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
1876  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
1877 
1878 #ifdef SELDON_LAPACK_CHECK_INFO
1879  if (info.GetInfo() != 0)
1880  throw LapackError(info.GetInfo(), "GetEigenvalues",
1881  "Failed to find eigenvalues ");
1882 #endif
1883 
1884  }
1885 
1886 
1887  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
1888  class Allocator3, class Allocator4>
1889  void GetEigenvaluesEigenvectors(Matrix<double, Prop1, RowSym, Allocator1>& A,
1890  Matrix<double, Prop2, RowSym, Allocator2>& B,
1891  Vector<double, VectFull, Allocator3>& w,
1892  Matrix<double, General, RowMajor, Allocator4>& z,
1893  LapackInfo& info)
1894  {
1895 #ifdef SELDON_CHECK_DIMENSIONS
1896  if (A.GetM() != B.GetM())
1897  throw WrongDim("GetEigenvaluesEigenvectors",
1898  "Matrix A and B must have the same size");
1899 #endif
1900 
1901  int itype = 1;
1902  int n = A.GetM();
1903  char uplo('L'), job('V');
1904  int lwork = 3*n;
1905  z.Reallocate(n,n);
1906  for (int i = 0; i < n; i++)
1907  for (int j = 0; j < n; j++)
1908  z(i,j) = A(i,j);
1909 
1910  A.Clear();
1911  Vector<double> work(lwork);
1912  w.Reallocate(n);
1913 
1914  dsygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
1915  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
1916 
1917 #ifdef SELDON_LAPACK_CHECK_INFO
1918  if (info.GetInfo() != 0)
1919  throw LapackError(info.GetInfo(), "GetEigenvalues",
1920  "Failed to find eigenvalues ");
1921 #endif
1922 
1923  Transpose(z);
1924  }
1925 
1926 
1927  template<class Prop1, class Prop2, class Allocator1,
1928  class Allocator2, class Allocator4, class Allocator5>
1929  void GetEigenvalues(Matrix<complex<double>, Prop1, RowSym, Allocator1>& A,
1930  Matrix<complex<double>, Prop2, RowSym, Allocator2>& B,
1931  Vector<complex<double>, VectFull, Allocator4>& alpha,
1932  Vector<complex<double>, VectFull, Allocator5>& beta,
1933  LapackInfo& info)
1934  {
1935 #ifdef SELDON_CHECK_DIMENSIONS
1936  if (A.GetM() != B.GetM())
1937  throw WrongDim("GetEigenvalues",
1938  "Matrix A and B must have the same size");
1939 #endif
1940 
1941  int n = A.GetM();
1942  Matrix<complex<double>, General, RowMajor, Allocator1> A2(n, n), B2(n, n);
1943  for (int i = 0; i < n; i++)
1944  for (int j = 0; j < n; j++)
1945  {
1946  A2(i, j) = A(i, j);
1947  B2(i, j) = B(i, j);
1948  }
1949 
1950  A.Clear();
1951  B.Clear();
1952  GetEigenvalues(A2, B2, alpha, beta, info);
1953  }
1954 
1955 
1956  template<class Prop1, class Prop2, class Prop3, class Allocator1,
1957  class Allocator2, class Allocator4,
1958  class Allocator5, class Allocator6>
1959  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
1960  Prop1, RowSym, Allocator1>& A,
1961  Matrix<complex<double>,
1962  Prop2, RowSym, Allocator2>& B,
1963  Vector<complex<double>,
1964  VectFull, Allocator4>& alpha,
1965  Vector<complex<double>,
1966  VectFull, Allocator5>& beta,
1967  Matrix<complex<double>,
1968  Prop3, RowMajor, Allocator6>& V,
1969  LapackInfo& info)
1970  {
1971 #ifdef SELDON_CHECK_DIMENSIONS
1972  if (A.GetM() != B.GetM())
1973  throw WrongDim("GetEigenvaluesEigenvectors",
1974  "Matrix A and B must have the same size");
1975 #endif
1976 
1977  int n = A.GetM();
1978  Matrix<complex<double>, General, RowMajor, Allocator1> A2(n, n), B2(n, n);
1979  for (int i = 0; i < n; i++)
1980  for (int j = 0; j < n; j++)
1981  {
1982  A2(i, j) = A(i, j);
1983  B2(i, j) = B(i, j);
1984  }
1985 
1986  A.Clear();
1987  B.Clear();
1988  GetEigenvaluesEigenvectors(A2, B2, alpha, beta, V, info);
1989  }
1990 
1991 
1992  /* ColSym */
1993 
1994 
1995  template<class Prop1, class Prop2, class Allocator1,
1996  class Allocator2, class Allocator3>
1997  void GetEigenvalues(Matrix<float, Prop1, ColSym, Allocator1>& A,
1998  Matrix<float, Prop2, ColSym, Allocator2>& B,
1999  Vector<float, VectFull, Allocator3>& w,
2000  LapackInfo& info)
2001  {
2002 #ifdef SELDON_CHECK_DIMENSIONS
2003  if (A.GetM() != B.GetM())
2004  throw WrongDim("GetEigenvalues",
2005  "Matrix A and B must have the same size");
2006 #endif
2007 
2008  int itype = 1;
2009  int n = A.GetM();
2010  char uplo('U'), job('N');
2011  int lwork = 3*n;
2012  Vector<float> work(lwork);
2013  w.Reallocate(n);
2014  ssygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2015  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
2016 
2017 #ifdef SELDON_LAPACK_CHECK_INFO
2018  if (info.GetInfo() != 0)
2019  throw LapackError(info.GetInfo(), "GetEigenvalues",
2020  "Failed to find eigenvalues ");
2021 #endif
2022 
2023  }
2024 
2025 
2026  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2027  class Allocator3, class Allocator4>
2028  void GetEigenvaluesEigenvectors(Matrix<float, Prop1, ColSym, Allocator1>& A,
2029  Matrix<float, Prop2, ColSym, Allocator2>& B,
2030  Vector<float, VectFull, Allocator3>& w,
2031  Matrix<float, General, ColMajor, Allocator4>& z,
2032  LapackInfo& info)
2033  {
2034 #ifdef SELDON_CHECK_DIMENSIONS
2035  if (A.GetM() != B.GetM())
2036  throw WrongDim("GetEigenvaluesEigenvectors",
2037  "Matrix A and B must have the same size");
2038 #endif
2039 
2040  int itype = 1;
2041  int n = A.GetM();
2042  char uplo('U'), job('V');
2043  int lwork = 3*n;
2044  Vector<float> work(lwork);
2045  w.Reallocate(n);
2046  z.Reallocate(n, n);
2047  for (int i = 0; i < n; i++)
2048  for (int j = 0; j < n; j++)
2049  z(i,j) = A(i,j);
2050 
2051  ssygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2052  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
2053 
2054 #ifdef SELDON_LAPACK_CHECK_INFO
2055  if (info.GetInfo() != 0)
2056  throw LapackError(info.GetInfo(), "GetEigenvalues",
2057  "Failed to find eigenvalues ");
2058 #endif
2059 
2060  }
2061 
2062 
2063  template<class Prop1, class Prop2, class Allocator1,
2064  class Allocator2, class Allocator4, class Allocator5>
2065  void GetEigenvalues(Matrix<complex<float>, Prop1, ColSym, Allocator1>& A,
2066  Matrix<complex<float>, Prop2, ColSym, Allocator2>& B,
2067  Vector<complex<float>, VectFull, Allocator4>& alpha,
2068  Vector<complex<float>, VectFull, Allocator5>& beta,
2069  LapackInfo& info)
2070  {
2071 #ifdef SELDON_CHECK_DIMENSIONS
2072  if (A.GetM() != B.GetM())
2073  throw WrongDim("GetEigenvalues",
2074  "Matrix A and B must have the same size");
2075 #endif
2076 
2077  int n = A.GetM();
2078  Matrix<complex<float>, General, ColMajor, Allocator1> A2(n, n), B2(n, n);
2079  for (int i = 0; i < n; i++)
2080  for (int j = 0; j < n; j++)
2081  {
2082  A2(i, j) = A(i, j);
2083  B2(i, j) = B(i, j);
2084  }
2085 
2086  A.Clear();
2087  B.Clear();
2088  GetEigenvalues(A2, B2, alpha, beta, info);
2089  }
2090 
2091 
2092  template<class Prop1, class Prop2, class Prop3, class Alloc1,
2093  class Alloc2, class Alloc4, class Alloc5, class Alloc6>
2094  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
2095  Prop1, ColSym, Alloc1>& A,
2096  Matrix<complex<float>,
2097  Prop2, ColSym, Alloc2>& B,
2098  Vector<complex<float>, VectFull, Alloc4>& alpha,
2099  Vector<complex<float>, VectFull, Alloc5>& beta,
2100  Matrix<complex<float>,
2101  Prop3, ColMajor, Alloc6>& V,
2102  LapackInfo& info)
2103  {
2104 #ifdef SELDON_CHECK_DIMENSIONS
2105  if (A.GetM() != B.GetM())
2106  throw WrongDim("GetEigenvaluesEigenvectors",
2107  "Matrix A and B must have the same size");
2108 #endif
2109 
2110  int n = A.GetM();
2111  Matrix<complex<float>, General, ColMajor, Alloc1> A2(n, n), B2(n, n);
2112  for (int i = 0; i < n; i++)
2113  for (int j = 0; j < n; j++)
2114  {
2115  A2(i, j) = A(i, j);
2116  B2(i, j) = B(i, j);
2117  }
2118 
2119  A.Clear();
2120  B.Clear();
2121  GetEigenvaluesEigenvectors(A2, B2, alpha, beta, V, info);
2122  }
2123 
2124  template<class Prop1, class Prop2, class Allocator1,
2125  class Allocator2, class Allocator3>
2126  void GetEigenvalues(Matrix<double, Prop1, ColSym, Allocator1>& A,
2127  Matrix<double, Prop2, ColSym, Allocator2>& B,
2128  Vector<double, VectFull, Allocator3>& w,
2129  LapackInfo& info)
2130  {
2131 #ifdef SELDON_CHECK_DIMENSIONS
2132  if (A.GetM() != B.GetM())
2133  throw WrongDim("GetEigenvalues",
2134  "Matrix A and B must have the same size");
2135 #endif
2136 
2137  int itype = 1;
2138  int n = A.GetM();
2139  char uplo('U'), job('N');
2140  int lwork = 3*n;
2141  Vector<double> work(lwork);
2142  w.Reallocate(n);
2143  dsygv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2144  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
2145 
2146 #ifdef SELDON_LAPACK_CHECK_INFO
2147  if (info.GetInfo() != 0)
2148  throw LapackError(info.GetInfo(), "GetEigenvalues",
2149  "Failed to find eigenvalues ");
2150 #endif
2151 
2152  }
2153 
2154 
2155  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2156  class Allocator3, class Allocator4>
2157  void GetEigenvaluesEigenvectors(Matrix<double, Prop1, ColSym, Allocator1>& A,
2158  Matrix<double, Prop2, ColSym, Allocator2>& B,
2159  Vector<double, VectFull, Allocator3>& w,
2160  Matrix<double, General, ColMajor, Allocator4>& z,
2161  LapackInfo& info)
2162  {
2163 #ifdef SELDON_CHECK_DIMENSIONS
2164  if (A.GetM() != B.GetM())
2165  throw WrongDim("GetEigenvaluesEigenvectors",
2166  "Matrix A and B must have the same size");
2167 #endif
2168 
2169  int itype = 1;
2170  int n = A.GetM();
2171  char uplo('U'), job('V');
2172  z.Reallocate(n,n);
2173  for (int i = 0; i < n; i++)
2174  for (int j = 0; j < n; j++)
2175  z(i,j) = A(i,j);
2176 
2177  A.Clear();
2178  int lwork = 3*n;
2179  Vector<double> work(lwork);
2180  w.Reallocate(n);
2181 
2182  dsygv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2183  w.GetData(), work.GetData(), &lwork, &info.GetInfoRef());
2184 
2185 #ifdef SELDON_LAPACK_CHECK_INFO
2186  if (info.GetInfo() != 0)
2187  throw LapackError(info.GetInfo(), "GetEigenvalues",
2188  "Failed to find eigenvalues ");
2189 #endif
2190 
2191  }
2192 
2193 
2194  template<class Prop1, class Prop2, class Allocator1,
2195  class Allocator2, class Allocator4, class Allocator5>
2196  void GetEigenvalues(Matrix<complex<double>, Prop1, ColSym, Allocator1>& A,
2197  Matrix<complex<double>, Prop2, ColSym, Allocator2>& B,
2198  Vector<complex<double>, VectFull, Allocator4>& alpha,
2199  Vector<complex<double>, VectFull, Allocator5>& beta,
2200  LapackInfo& info)
2201  {
2202 #ifdef SELDON_CHECK_DIMENSIONS
2203  if (A.GetM() != B.GetM())
2204  throw WrongDim("GetEigenvalues",
2205  "Matrix A and B must have the same size");
2206 #endif
2207 
2208  int n = A.GetM();
2209  Matrix<complex<double>, General, ColMajor, Allocator1> A2(n, n), B2(n, n);
2210  for (int i = 0; i < n; i++)
2211  for (int j = 0; j < n; j++)
2212  {
2213  A2(i, j) = A(i, j);
2214  B2(i, j) = B(i, j);
2215  }
2216 
2217  A.Clear();
2218  B.Clear();
2219  GetEigenvalues(A2, B2, alpha, beta, info);
2220  }
2221 
2222 
2223  template<class Prop1, class Prop2, class Prop3, class Alloc1,
2224  class Alloc2, class Alloc4, class Alloc5, class Alloc6>
2225  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
2226  Prop1, ColSym, Alloc1>& A,
2227  Matrix<complex<double>,
2228  Prop2, ColSym, Alloc2>& B,
2229  Vector<complex<double>, VectFull, Alloc4>& alpha,
2230  Vector<complex<double>, VectFull, Alloc5>& beta,
2231  Matrix<complex<double>,
2232  Prop3, ColMajor, Alloc6>& V,
2233  LapackInfo& info)
2234  {
2235 #ifdef SELDON_CHECK_DIMENSIONS
2236  if (A.GetM() != B.GetM())
2237  throw WrongDim("GetEigenvaluesEigenvectors",
2238  "Matrix A and B must have the same size");
2239 #endif
2240 
2241  int n = A.GetM();
2242  Matrix<complex<double>, General, ColMajor, Alloc1> A2(n, n), B2(n, n);
2243  for (int i = 0; i < n; i++)
2244  for (int j = 0; j < n; j++)
2245  {
2246  A2(i, j) = A(i, j);
2247  B2(i, j) = B(i, j);
2248  }
2249 
2250  A.Clear();
2251  B.Clear();
2252  GetEigenvaluesEigenvectors(A2, B2, alpha, beta, V, info);
2253  }
2254 
2255 
2256  /* RowHerm */
2257 
2258 
2259  template<class Prop1, class Prop2, class Allocator1,
2260  class Allocator2, class Allocator3>
2261  void GetEigenvalues(Matrix<complex<float>, Prop1, RowHerm, Allocator1>& A,
2262  Matrix<complex<float>, Prop2, RowHerm, Allocator2>& B,
2263  Vector<float, VectFull, Allocator3>& w,
2264  LapackInfo& info)
2265  {
2266 #ifdef SELDON_CHECK_DIMENSIONS
2267  if (A.GetM() != B.GetM())
2268  throw WrongDim("GetEigenvalues",
2269  "Matrix A and B must have the same size");
2270 #endif
2271 
2272  int itype = 1;
2273  int n = A.GetM();
2274  char uplo('L'), job('N');
2275  int lwork = 2*n;
2276  Vector<complex<float> > work(lwork);
2277  Vector<float> rwork(3*n);
2278  w.Reallocate(n);
2279  Conjugate(A); Conjugate(B);
2280  chegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2281  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2282  &info.GetInfoRef());
2283 
2284 #ifdef SELDON_LAPACK_CHECK_INFO
2285  if (info.GetInfo() != 0)
2286  throw LapackError(info.GetInfo(), "GetEigenvalues",
2287  "Failed to find eigenvalues ");
2288 #endif
2289 
2290  }
2291 
2292 
2293  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2294  class Allocator3, class Allocator4>
2295  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
2296  Prop1, RowHerm, Allocator1>& A,
2297  Matrix<complex<float>,
2298  Prop2, RowHerm, Allocator2>& B,
2299  Vector<float, VectFull, Allocator3>& w,
2300  Matrix<complex<float>,
2301  General, RowMajor, Allocator4>& z,
2302  LapackInfo& info)
2303  {
2304 #ifdef SELDON_CHECK_DIMENSIONS
2305  if (A.GetM() != B.GetM())
2306  throw WrongDim("GetEigenvaluesEigenvectors",
2307  "Matrix A and B must have the same size");
2308 #endif
2309 
2310  int itype = 1;
2311  int n = A.GetM();
2312  char uplo('L'), job('V');
2313  int lwork = 2*n;
2314  z.Reallocate(n,n);
2315  for (int i = 0; i < n; i++)
2316  for (int j = 0; j < n; j++)
2317  z(i,j) = conj(A(i,j));
2318 
2319  A.Clear();
2320  Vector<complex<float> > work(lwork);
2321  w.Reallocate(n);
2322  Vector<float> rwork(3*n);
2323  Conjugate(B);
2324 
2325  chegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2326  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2327  &info.GetInfoRef());
2328 
2329 #ifdef SELDON_LAPACK_CHECK_INFO
2330  if (info.GetInfo() != 0)
2331  throw LapackError(info.GetInfo(), "GetEigenvalues",
2332  "Failed to find eigenvalues ");
2333 #endif
2334 
2335  Transpose(z);
2336  }
2337 
2338 
2339  template<class Prop1, class Prop2, class Allocator1,
2340  class Allocator2, class Allocator3>
2341  void GetEigenvalues(Matrix<complex<double>, Prop1, RowHerm, Allocator1>& A,
2342  Matrix<complex<double>, Prop2, RowHerm, Allocator2>& B,
2343  Vector<double, VectFull, Allocator3>& w,
2344  LapackInfo& info)
2345  {
2346 #ifdef SELDON_CHECK_DIMENSIONS
2347  if (A.GetM() != B.GetM())
2348  throw WrongDim("GetEigenvalues",
2349  "Matrix A and B must have the same size");
2350 #endif
2351 
2352  int itype = 1;
2353  int n = A.GetM();
2354  char uplo('L'), job('N');
2355 
2356  int lwork = 3*n;
2357  Vector<complex<double> > work(lwork);
2358  Vector<double> rwork(3*n);
2359  w.Reallocate(n);
2360  Conjugate(A);
2361  Conjugate(B);
2362 
2363  zhegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2364  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2365  &info.GetInfoRef());
2366 
2367 #ifdef SELDON_LAPACK_CHECK_INFO
2368  if (info.GetInfo() != 0)
2369  throw LapackError(info.GetInfo(), "GetEigenvalues",
2370  "Failed to find eigenvalues ");
2371 #endif
2372 
2373  }
2374 
2375 
2376  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2377  class Allocator3, class Allocator4>
2378  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
2379  Prop1, RowHerm, Allocator1>& A,
2380  Matrix<complex<double>,
2381  Prop2, RowHerm, Allocator2>& B,
2382  Vector<double, VectFull, Allocator3>& w,
2383  Matrix<complex<double>,
2384  General, RowMajor, Allocator4>& z,
2385  LapackInfo& info)
2386  {
2387 #ifdef SELDON_CHECK_DIMENSIONS
2388  if (A.GetM() != B.GetM())
2389  throw WrongDim("GetEigenvaluesEigenvectors",
2390  "Matrix A and B must have the same size");
2391 #endif
2392 
2393  int itype = 1;
2394  int n = A.GetM();
2395  char uplo('L'), job('V');
2396  z.Reallocate(n,n);
2397  for (int i = 0; i < n; i++)
2398  for (int j = 0; j < n; j++)
2399  z(i,j) = conj(A(i,j));
2400 
2401  A.Clear();
2402  int lwork = 3*n;
2403  Vector<complex<double> > work(lwork);
2404  w.Reallocate(n);
2405  Conjugate(B);
2406 
2407  Vector<double> rwork(3*n);
2408  zhegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2409  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2410  &info.GetInfoRef());
2411 
2412 #ifdef SELDON_LAPACK_CHECK_INFO
2413  if (info.GetInfo() != 0)
2414  throw LapackError(info.GetInfo(), "GetEigenvalues",
2415  "Failed to find eigenvalues ");
2416 #endif
2417 
2418  Transpose(z);
2419  }
2420 
2421 
2422  /* ColHerm */
2423 
2424 
2425  template<class Prop1, class Prop2, class Allocator1,
2426  class Allocator2, class Allocator3>
2427  void GetEigenvalues(Matrix<complex<float>, Prop1, ColHerm, Allocator1>& A,
2428  Matrix<complex<float>, Prop2, ColHerm, Allocator2>& B,
2429  Vector<float, VectFull, Allocator3>& w,
2430  LapackInfo& info)
2431  {
2432 #ifdef SELDON_CHECK_DIMENSIONS
2433  if (A.GetM() != B.GetM())
2434  throw WrongDim("GetEigenvalues",
2435  "Matrix A and B must have the same size");
2436 #endif
2437 
2438  int itype = 1;
2439  int n = A.GetM();
2440  char uplo('U'), job('N');
2441  int lwork = 2*n;
2442  Vector<complex<float> > work(lwork);
2443  Vector<float> rwork(3*n);
2444  w.Reallocate(n);
2445  chegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2446  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2447  &info.GetInfoRef());
2448 
2449 #ifdef SELDON_LAPACK_CHECK_INFO
2450  if (info.GetInfo() != 0)
2451  throw LapackError(info.GetInfo(), "GetEigenvalues",
2452  "Failed to find eigenvalues ");
2453 #endif
2454 
2455  }
2456 
2457 
2458  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2459  class Allocator3, class Allocator4>
2460  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
2461  Prop1, ColHerm, Allocator1>& A,
2462  Matrix<complex<float>,
2463  Prop2, ColHerm, Allocator2>& B,
2464  Vector<float, VectFull, Allocator3>& w,
2465  Matrix<complex<float>,
2466  General, ColMajor, Allocator4>& z,
2467  LapackInfo& info)
2468  {
2469 #ifdef SELDON_CHECK_DIMENSIONS
2470  if (A.GetM() != B.GetM())
2471  throw WrongDim("GetEigenvaluesEigenvectors",
2472  "Matrix A and B must have the same size");
2473 #endif
2474 
2475  int itype = 1;
2476  int n = A.GetM();
2477  char uplo('U'); char job('V');
2478  z.Reallocate(n,n);
2479  for (int i = 0; i < n; i++)
2480  for (int j = 0; j < n; j++)
2481  z(i,j) = A(i,j);
2482 
2483  A.Clear();
2484  int lwork = 3*n;
2485  Vector<complex<float> > work(lwork);
2486  w.Reallocate(n);
2487  Vector<float> rwork(3*n);
2488  chegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2489  w.GetData(), work.GetData(), &lwork,
2490  rwork.GetData(), &info.GetInfoRef());
2491 
2492 #ifdef SELDON_LAPACK_CHECK_INFO
2493  if (info.GetInfo() != 0)
2494  throw LapackError(info.GetInfo(), "GetEigenvalues",
2495  "Failed to find eigenvalues ");
2496 #endif
2497 
2498  }
2499 
2500 
2501  template<class Prop1, class Prop2, class Allocator1,
2502  class Allocator2, class Allocator3>
2503  void GetEigenvalues(Matrix<complex<double>, Prop1, ColHerm, Allocator1>& A,
2504  Matrix<complex<double>, Prop2, ColHerm, Allocator2>& B,
2505  Vector<double, VectFull, Allocator3>& w,
2506  LapackInfo& info)
2507  {
2508 #ifdef SELDON_CHECK_DIMENSIONS
2509  if (A.GetM() != B.GetM())
2510  throw WrongDim("GetEigenvalues",
2511  "Matrix A and B must have the same size");
2512 #endif
2513 
2514  int itype = 1;
2515  int n = A.GetM();
2516  char uplo('U'), job('N');
2517  int lwork = 3*n;
2518  Vector<complex<double> > work(lwork);
2519  Vector<double> rwork(3*n);
2520  w.Reallocate(n);
2521  zhegv_(&itype, &job, &uplo, &n, A.GetData(), &n, B.GetData(), &n,
2522  w.GetData(), work.GetData(), &lwork, rwork.GetData(),
2523  &info.GetInfoRef());
2524 
2525 #ifdef SELDON_LAPACK_CHECK_INFO
2526  if (info.GetInfo() != 0)
2527  throw LapackError(info.GetInfo(), "GetEigenvalues",
2528  "Failed to find eigenvalues ");
2529 #endif
2530 
2531  }
2532 
2533 
2534  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2535  class Allocator3, class Allocator4>
2536  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
2537  Prop1, ColHerm, Allocator1>& A,
2538  Matrix<complex<double>,
2539  Prop2, ColHerm, Allocator2>& B,
2540  Vector<double, VectFull, Allocator3>& w,
2541  Matrix<complex<double>,
2542  General, ColMajor, Allocator4>& z,
2543  LapackInfo& info)
2544  {
2545 #ifdef SELDON_CHECK_DIMENSIONS
2546  if (A.GetM() != B.GetM())
2547  throw WrongDim("GetEigenvaluesEigenvectors",
2548  "Matrix A and B must have the same size");
2549 #endif
2550 
2551  int itype = 1;
2552  int n = A.GetM();
2553  char uplo('U'), job('V');
2554  z.Reallocate(n,n);
2555  for (int i = 0; i < n; i++)
2556  for (int j = 0; j < n; j++)
2557  z(i,j) = A(i,j);
2558 
2559  A.Clear();
2560  int lwork = 3*n;
2561  Vector<complex<double> > work(lwork);
2562  w.Reallocate(n);
2563  Vector<double> rwork(3*n);
2564 
2565  zhegv_(&itype, &job, &uplo, &n, z.GetData(), &n, B.GetData(), &n,
2566  w.GetData(), work.GetData(), &lwork,
2567  rwork.GetData(), &info.GetInfoRef());
2568 
2569 #ifdef SELDON_LAPACK_CHECK_INFO
2570  if (info.GetInfo() != 0)
2571  throw LapackError(info.GetInfo(), "GetEigenvalues",
2572  "Failed to find eigenvalues ");
2573 #endif
2574 
2575  }
2576 
2577 
2578  /* RowSymPacked */
2579 
2580 
2581  template<class Prop1, class Prop2, class Allocator1,
2582  class Allocator2, class Allocator3>
2583  void GetEigenvalues(Matrix<float, Prop1, RowSymPacked, Allocator1>& A,
2584  Matrix<float, Prop2, RowSymPacked, Allocator2>& B,
2585  Vector<float, VectFull, Allocator3>& w,
2586  LapackInfo& info)
2587  {
2588 #ifdef SELDON_CHECK_DIMENSIONS
2589  if (A.GetM() != B.GetM())
2590  throw WrongDim("GetEigenvalues",
2591  "Matrix A and B must have the same size");
2592 #endif
2593 
2594  int itype = 1;
2595  int n = A.GetM();
2596  char uplo('L'); char job('N');
2597  int lwork = 3*n; Vector<float> work(lwork);
2598  w.Reallocate(n);
2599  sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2600  A.GetData(), &n, work.GetData(), &info.GetInfoRef());
2601 
2602 #ifdef SELDON_LAPACK_CHECK_INFO
2603  if (info.GetInfo() != 0)
2604  throw LapackError(info.GetInfo(), "GetEigenvalues",
2605  "Failed to find eigenvalues ");
2606 #endif
2607 
2608  }
2609 
2610 
2611  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2612  class Allocator3, class Allocator4>
2613  void GetEigenvaluesEigenvectors(Matrix<float,
2614  Prop1, RowSymPacked, Allocator1>& A,
2615  Matrix<float,
2616  Prop2, RowSymPacked, Allocator2>& B,
2617  Vector<float, VectFull, Allocator3>& w,
2618  Matrix<float,
2619  General, RowMajor, Allocator4>& z,
2620  LapackInfo& info)
2621  {
2622 #ifdef SELDON_CHECK_DIMENSIONS
2623  if (A.GetM() != B.GetM())
2624  throw WrongDim("GetEigenvaluesEigenvectors",
2625  "Matrix A and B must have the same size");
2626 #endif
2627 
2628  int itype = 1;
2629  int n = A.GetM();
2630  char uplo('L'); char job('V');
2631  int lwork = 3*n; Vector<float> work(lwork);
2632  w.Reallocate(n);
2633  z.Reallocate(n,n);
2634  sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2635  z.GetData(), &n, work.GetData(), &info.GetInfoRef());
2636 
2637 #ifdef SELDON_LAPACK_CHECK_INFO
2638  if (info.GetInfo() != 0)
2639  throw LapackError(info.GetInfo(), "GetEigenvalues",
2640  "Failed to find eigenvalues ");
2641 #endif
2642 
2643  Transpose(z);
2644  }
2645 
2646 
2647  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2648  class Allocator3, class Allocator4>
2649  void GetEigenvalues(Matrix<complex<float>,
2650  Prop1, RowSymPacked, Allocator1>& A,
2651  Matrix<complex<float>,
2652  Prop2, RowSymPacked, Allocator2>& B,
2653  Vector<complex<float>, VectFull, Allocator3>& alpha,
2654  Vector<complex<float>, VectFull, Allocator4>& beta,
2655  LapackInfo& info)
2656  {
2657 #ifdef SELDON_CHECK_DIMENSIONS
2658  if (A.GetM() != B.GetM())
2659  throw WrongDim("GetEigenvalues",
2660  "Matrix A and B must have the same size");
2661 #endif
2662 
2663  int n = A.GetM();
2664  Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n);
2665  for (int i = 0; i < n; i++)
2666  for (int j = 0; j < n; j++)
2667  {
2668  C(i,j) = A(i,j);
2669  D(i,j) = B(i,j);
2670  }
2671 
2672  A.Clear();
2673  B.Clear();
2674  alpha.Reallocate(n);
2675  beta.Reallocate(n);
2676  GetEigenvalues(C, D, alpha, beta, info);
2677  }
2678 
2679 
2680  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2681  class Allocator3, class Allocator4, class Allocator5>
2682  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
2683  Prop1, RowSymPacked, Allocator1>& A,
2684  Matrix<complex<float>,
2685  Prop2, RowSymPacked, Allocator2>& B,
2686  Vector<complex<float>,
2687  VectFull, Allocator3>& alpha,
2688  Vector<complex<float>,
2689  VectFull, Allocator4>& beta,
2690  Matrix<complex<float>,
2691  General, RowMajor, Allocator5>& z,
2692  LapackInfo& info)
2693  {
2694 #ifdef SELDON_CHECK_DIMENSIONS
2695  if (A.GetM() != B.GetM())
2696  throw WrongDim("GetEigenvaluesEigenvectors",
2697  "Matrix A and B must have the same size");
2698 #endif
2699 
2700  int n = A.GetM();
2701  Matrix<complex<float>, General, RowMajor> C(n,n), D(n,n);
2702  for (int i = 0; i < n; i++)
2703  for (int j = 0; j < n; j++)
2704  {
2705  C(i,j) = A(i,j);
2706  D(i,j) = B(i,j);
2707  }
2708 
2709  A.Clear();
2710  B.Clear();
2711  alpha.Reallocate(n);
2712  beta.Reallocate(n);
2713  z.Reallocate(n,n);
2714  GetEigenvaluesEigenvectors(C, D, alpha, beta, z, info);
2715  }
2716 
2717 
2718  template<class Prop1, class Prop2, class Allocator1,
2719  class Allocator2, class Allocator3>
2720  void GetEigenvalues(Matrix<double, Prop1, RowSymPacked, Allocator1>& A,
2721  Matrix<double, Prop2, RowSymPacked, Allocator2>& B,
2722  Vector<double, VectFull, Allocator3>& w,
2723  LapackInfo& info)
2724  {
2725 #ifdef SELDON_CHECK_DIMENSIONS
2726  if (A.GetM() != B.GetM())
2727  throw WrongDim("GetEigenvalues",
2728  "Matrix A and B must have the same size");
2729 #endif
2730 
2731  int itype = 1;
2732  int n = A.GetM();
2733  char uplo('L'); char job('N');
2734  int lwork = 3*n; Vector<double> work(lwork);
2735  w.Reallocate(n);
2736  dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2737  A.GetData(), &n, work.GetData(), &info.GetInfoRef());
2738 
2739 #ifdef SELDON_LAPACK_CHECK_INFO
2740  if (info.GetInfo() != 0)
2741  throw LapackError(info.GetInfo(), "GetEigenvalues",
2742  "Failed to find eigenvalues ");
2743 #endif
2744 
2745  }
2746 
2747 
2748  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2749  class Allocator3, class Allocator4>
2750  void GetEigenvaluesEigenvectors(Matrix<double,
2751  Prop1, RowSymPacked, Allocator1>& A,
2752  Matrix<double,
2753  Prop2, RowSymPacked, Allocator2>& B,
2754  Vector<double, VectFull, Allocator3>& w,
2755  Matrix<double,
2756  General, RowMajor, Allocator4>& z,
2757  LapackInfo& info)
2758  {
2759 #ifdef SELDON_CHECK_DIMENSIONS
2760  if (A.GetM() != B.GetM())
2761  throw WrongDim("GetEigenvaluesEigenvectors",
2762  "Matrix A and B must have the same size");
2763 #endif
2764 
2765  int itype = 1;
2766  int n = A.GetM();
2767  char uplo('L'); char job('V');
2768  int lwork = 3*n; Vector<double> work(lwork);
2769  w.Reallocate(n);
2770  z.Reallocate(n,n);
2771  dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2772  z.GetData(), &n, work.GetData(), &info.GetInfoRef());
2773 
2774 #ifdef SELDON_LAPACK_CHECK_INFO
2775  if (info.GetInfo() != 0)
2776  throw LapackError(info.GetInfo(), "GetEigenvalues",
2777  "Failed to find eigenvalues ");
2778 #endif
2779 
2780  Transpose(z);
2781  }
2782 
2783 
2784  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2785  class Allocator3, class Allocator4>
2786  void GetEigenvalues(Matrix<complex<double>,
2787  Prop1, RowSymPacked, Allocator1>& A,
2788  Matrix<complex<double>,
2789  Prop2, RowSymPacked, Allocator2>& B,
2790  Vector<complex<double>, VectFull, Allocator3>& alpha,
2791  Vector<complex<double>, VectFull, Allocator4>& beta,
2792  LapackInfo& info)
2793  {
2794 #ifdef SELDON_CHECK_DIMENSIONS
2795  if (A.GetM() != B.GetM())
2796  throw WrongDim("GetEigenvalues",
2797  "Matrix A and B must have the same size");
2798 #endif
2799 
2800  int n = A.GetM();
2801  Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n);
2802  for (int i = 0; i < n; i++)
2803  for (int j = 0; j < n; j++)
2804  {
2805  C(i,j) = A(i,j);
2806  D(i,j) = B(i,j);
2807  }
2808 
2809  A.Clear();
2810  B.Clear();
2811  alpha.Reallocate(n);
2812  beta.Reallocate(n);
2813  GetEigenvalues(C, D, alpha, beta, info);
2814  }
2815 
2816 
2817  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2818  class Allocator3, class Allocator4, class Allocator5>
2819  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
2820  Prop1, RowSymPacked, Allocator1>& A,
2821  Matrix<complex<double>,
2822  Prop2, RowSymPacked, Allocator2>& B,
2823  Vector<complex<double>,
2824  VectFull, Allocator3>& alpha,
2825  Vector<complex<double>,
2826  VectFull, Allocator4>& beta,
2827  Matrix<complex<double>,
2828  General, RowMajor, Allocator5>& z,
2829  LapackInfo& info)
2830  {
2831 #ifdef SELDON_CHECK_DIMENSIONS
2832  if (A.GetM() != B.GetM())
2833  throw WrongDim("GetEigenvaluesEigenvectors",
2834  "Matrix A and B must have the same size");
2835 #endif
2836 
2837  int n = A.GetM();
2838  Matrix<complex<double>, General, RowMajor> C(n,n), D(n,n);
2839  for (int i = 0; i < n; i++)
2840  for (int j = 0; j < n; j++)
2841  {
2842  C(i,j) = A(i,j);
2843  D(i,j) = B(i,j);
2844  }
2845 
2846  A.Clear();
2847  B.Clear();
2848  alpha.Reallocate(n);
2849  beta.Reallocate(n);
2850  z.Reallocate(n,n);
2851  GetEigenvaluesEigenvectors(C, D, alpha, beta, z, info);
2852  }
2853 
2854 
2855  /* ColSymPacked */
2856 
2857 
2858  template<class Prop1, class Prop2, class Allocator1,
2859  class Allocator2, class Allocator3>
2860  void GetEigenvalues(Matrix<float, Prop1, ColSymPacked, Allocator1>& A,
2861  Matrix<float, Prop2, ColSymPacked, Allocator2>& B,
2862  Vector<float, VectFull, Allocator3>& w,
2863  LapackInfo& info)
2864  {
2865 #ifdef SELDON_CHECK_DIMENSIONS
2866  if (A.GetM() != B.GetM())
2867  throw WrongDim("GetEigenvalues",
2868  "Matrix A and B must have the same size");
2869 #endif
2870 
2871  int itype = 1;
2872  int n = A.GetM();
2873  char uplo('U'); char job('N');
2874  int lwork = 3*n; Vector<float> work(lwork);
2875  w.Reallocate(n);
2876  sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2877  A.GetData(), &n, work.GetData(), &info.GetInfoRef());
2878 
2879 #ifdef SELDON_LAPACK_CHECK_INFO
2880  if (info.GetInfo() != 0)
2881  throw LapackError(info.GetInfo(), "GetEigenvalues",
2882  "Failed to find eigenvalues ");
2883 #endif
2884 
2885  }
2886 
2887 
2888  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2889  class Allocator3, class Allocator4>
2890  void GetEigenvaluesEigenvectors(Matrix<float,
2891  Prop1, ColSymPacked, Allocator1>& A,
2892  Matrix<float,
2893  Prop2, ColSymPacked, Allocator2>& B,
2894  Vector<float, VectFull, Allocator3>& w,
2895  Matrix<float,
2896  General, ColMajor, Allocator4>& z,
2897  LapackInfo& info)
2898  {
2899 #ifdef SELDON_CHECK_DIMENSIONS
2900  if (A.GetM() != B.GetM())
2901  throw WrongDim("GetEigenvaluesEigenvectors",
2902  "Matrix A and B must have the same size");
2903 #endif
2904 
2905  int itype = 1;
2906  int n = A.GetM();
2907  char uplo('U'); char job('V'); int lwork = 3*n;
2908  Vector<float> work(lwork);
2909  w.Reallocate(n);
2910  z.Reallocate(n,n);
2911  sspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
2912  z.GetData(), &n, work.GetData(), &info.GetInfoRef());
2913 
2914 #ifdef SELDON_LAPACK_CHECK_INFO
2915  if (info.GetInfo() != 0)
2916  throw LapackError(info.GetInfo(), "GetEigenvalues",
2917  "Failed to find eigenvalues ");
2918 #endif
2919 
2920  }
2921 
2922 
2923  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2924  class Allocator3, class Allocator4>
2925  void GetEigenvalues(Matrix<complex<float>,
2926  Prop1, ColSymPacked, Allocator1>& A,
2927  Matrix<complex<float>,
2928  Prop2, ColSymPacked, Allocator2>& B,
2929  Vector<complex<float>, VectFull, Allocator3>& alpha,
2930  Vector<complex<float>, VectFull, Allocator4>& beta,
2931  LapackInfo& info)
2932  {
2933 #ifdef SELDON_CHECK_DIMENSIONS
2934  if (A.GetM() != B.GetM())
2935  throw WrongDim("GetEigenvalues",
2936  "Matrix A and B must have the same size");
2937 #endif
2938 
2939  int n = A.GetM();
2940  Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n);
2941  for (int i = 0; i < n; i++)
2942  for (int j = 0; j < n; j++)
2943  {
2944  C(i,j) = A(i,j);
2945  D(i,j) = B(i,j);
2946  }
2947 
2948  A.Clear();
2949  B.Clear();
2950  alpha.Reallocate(n);
2951  beta.Reallocate(n);
2952  GetEigenvalues(C, D, alpha, beta, info);
2953  }
2954 
2955 
2956  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
2957  class Allocator3, class Allocator4, class Allocator5>
2958  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
2959  Prop1, ColSymPacked, Allocator1>& A,
2960  Matrix<complex<float>,
2961  Prop2, ColSymPacked, Allocator2>& B,
2962  Vector<complex<float>,
2963  VectFull, Allocator3>& alpha,
2964  Vector<complex<float>,
2965  VectFull, Allocator4>& beta,
2966  Matrix<complex<float>,
2967  General, ColMajor, Allocator5>& z,
2968  LapackInfo& info)
2969  {
2970 #ifdef SELDON_CHECK_DIMENSIONS
2971  if (A.GetM() != B.GetM())
2972  throw WrongDim("GetEigenvaluesEigenvectors",
2973  "Matrix A and B must have the same size");
2974 #endif
2975 
2976  int n = A.GetM();
2977  Matrix<complex<float>, General, ColMajor> C(n,n), D(n,n);
2978  for (int i = 0; i < n; i++)
2979  for (int j = 0; j < n; j++)
2980  {
2981  C(i,j) = A(i,j);
2982  D(i,j) = B(i,j);
2983  }
2984 
2985  A.Clear();
2986  B.Clear();
2987  alpha.Reallocate(n);
2988  beta.Reallocate(n);
2989  z.Reallocate(n,n);
2990  GetEigenvaluesEigenvectors(C, D, alpha, beta, z, info);
2991  }
2992 
2993 
2994  template<class Prop1, class Prop2, class Allocator1,
2995  class Allocator2, class Allocator3>
2996  void GetEigenvalues(Matrix<double, Prop1, ColSymPacked, Allocator1>& A,
2997  Matrix<double, Prop2, ColSymPacked, Allocator2>& B,
2998  Vector<double, VectFull, Allocator3>& w,
2999  LapackInfo& info)
3000  {
3001 #ifdef SELDON_CHECK_DIMENSIONS
3002  if (A.GetM() != B.GetM())
3003  throw WrongDim("GetEigenvalues",
3004  "Matrix A and B must have the same size");
3005 #endif
3006 
3007  int itype = 1;
3008  int n = A.GetM();
3009  char uplo('U'); char job('N');
3010  int lwork = 3*n; Vector<double> work(lwork);
3011  w.Reallocate(n);
3012  dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3013  A.GetData(), &n, work.GetData(), &info.GetInfoRef());
3014 
3015 #ifdef SELDON_LAPACK_CHECK_INFO
3016  if (info.GetInfo() != 0)
3017  throw LapackError(info.GetInfo(), "GetEigenvalues",
3018  "Failed to find eigenvalues ");
3019 #endif
3020 
3021  }
3022 
3023 
3024  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3025  class Allocator3, class Allocator4>
3026  void GetEigenvaluesEigenvectors(Matrix<double,
3027  Prop1, ColSymPacked, Allocator1>& A,
3028  Matrix<double,
3029  Prop2, ColSymPacked, Allocator2>& B,
3030  Vector<double, VectFull, Allocator3>& w,
3031  Matrix<double,
3032  General, ColMajor, Allocator4>& z,
3033  LapackInfo& info)
3034  {
3035 #ifdef SELDON_CHECK_DIMENSIONS
3036  if (A.GetM() != B.GetM())
3037  throw WrongDim("GetEigenvaluesEigenvectors",
3038  "Matrix A and B must have the same size");
3039 #endif
3040 
3041  int itype = 1;
3042  int n = A.GetM();
3043  char uplo('U'); char job('V'); int lwork = 3*n;
3044  Vector<double> work(lwork);
3045  w.Reallocate(n);
3046  z.Reallocate(n,n);
3047  dspgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3048  z.GetData(), &n, work.GetData(), &info.GetInfoRef());
3049 
3050 #ifdef SELDON_LAPACK_CHECK_INFO
3051  if (info.GetInfo() != 0)
3052  throw LapackError(info.GetInfo(), "GetEigenvalues",
3053  "Failed to find eigenvalues ");
3054 #endif
3055 
3056  }
3057 
3058 
3059  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3060  class Allocator3, class Allocator4>
3061  void GetEigenvalues(Matrix<complex<double>,
3062  Prop1, ColSymPacked, Allocator1>& A,
3063  Matrix<complex<double>,
3064  Prop2, ColSymPacked, Allocator2>& B,
3065  Vector<complex<double>, VectFull, Allocator3>& alpha,
3066  Vector<complex<double>, VectFull, Allocator4>& beta,
3067  LapackInfo& info)
3068  {
3069 #ifdef SELDON_CHECK_DIMENSIONS
3070  if (A.GetM() != B.GetM())
3071  throw WrongDim("GetEigenvalues",
3072  "Matrix A and B must have the same size");
3073 #endif
3074 
3075  int n = A.GetM();
3076  Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n);
3077  for (int i = 0; i < n; i++)
3078  for (int j = 0; j < n; j++)
3079  {
3080  C(i,j) = A(i,j);
3081  D(i,j) = B(i,j);
3082  }
3083 
3084  A.Clear();
3085  B.Clear();
3086  alpha.Reallocate(n);
3087  beta.Reallocate(n);
3088  GetEigenvalues(C, D, alpha, beta, info);
3089  }
3090 
3091 
3092  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3093  class Allocator3, class Allocator4, class Allocator5>
3094  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
3095  Prop1, ColSymPacked, Allocator1>& A,
3096  Matrix<complex<double>,
3097  Prop2, ColSymPacked, Allocator2>& B,
3098  Vector<complex<double>,
3099  VectFull, Allocator3>& alpha,
3100  Vector<complex<double>,
3101  VectFull, Allocator4>& beta,
3102  Matrix<complex<double>,
3103  General, ColMajor, Allocator5>& z,
3104  LapackInfo& info)
3105  {
3106 #ifdef SELDON_CHECK_DIMENSIONS
3107  if (A.GetM() != B.GetM())
3108  throw WrongDim("GetEigenvaluesEigenvectors",
3109  "Matrix A and B must have the same size");
3110 #endif
3111 
3112  int n = A.GetM();
3113  Matrix<complex<double>, General, ColMajor> C(n,n), D(n,n);
3114  for (int i = 0; i < n; i++)
3115  for (int j = 0; j < n; j++)
3116  {
3117  C(i,j) = A(i,j);
3118  D(i,j) = B(i,j);
3119  }
3120 
3121  A.Clear();
3122  B.Clear();
3123  alpha.Reallocate(n);
3124  beta.Reallocate(n);
3125  z.Reallocate(n,n);
3126  GetEigenvaluesEigenvectors(C, D, alpha, beta, z, info);
3127  }
3128 
3129 
3130  /* RowHermPacked */
3131 
3132 
3133  template<class Prop1, class Prop2, class Allocator1,
3134  class Allocator2, class Allocator3>
3135  void GetEigenvalues(Matrix<complex<float>,
3136  Prop1, RowHermPacked, Allocator1>& A,
3137  Matrix<complex<float>,
3138  Prop2, RowHermPacked, Allocator2>& B,
3139  Vector<float, VectFull, Allocator3>& w,
3140  LapackInfo& info)
3141  {
3142 #ifdef SELDON_CHECK_DIMENSIONS
3143  if (A.GetM() != B.GetM())
3144  throw WrongDim("GetEigenvalues",
3145  "Matrix A and B must have the same size");
3146 #endif
3147 
3148  int itype = 1;
3149  int n = A.GetM();
3150  char uplo('L'), job('N');
3151  int lwork = 2*n;
3152  Vector<complex<float> > work(lwork);
3153  Vector<float> rwork(3*n);
3154  w.Reallocate(n);
3155  Conjugate(A); Conjugate(B);
3156  chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3157  A.GetData(), &n, work.GetData(), rwork.GetData(),
3158  &info.GetInfoRef());
3159 
3160 #ifdef SELDON_LAPACK_CHECK_INFO
3161  if (info.GetInfo() != 0)
3162  throw LapackError(info.GetInfo(), "GetEigenvalues",
3163  "Failed to find eigenvalues ");
3164 #endif
3165 
3166  }
3167 
3168 
3169  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3170  class Allocator3, class Allocator4>
3171  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
3172  Prop1, RowHermPacked, Allocator1>& A,
3173  Matrix<complex<float>,
3174  Prop2, RowHermPacked, Allocator2>& B,
3175  Vector<float, VectFull, Allocator3>& w,
3176  Matrix<complex<float>,
3177  General, RowMajor, Allocator4>& z,
3178  LapackInfo& info)
3179  {
3180 #ifdef SELDON_CHECK_DIMENSIONS
3181  if (A.GetM() != B.GetM())
3182  throw WrongDim("GetEigenvaluesEigenvectors",
3183  "Matrix A and B must have the same size");
3184 #endif
3185 
3186  int itype = 1;
3187  int n = A.GetM();
3188  char uplo('L'); char job('V'); int lwork = 2*n;
3189  Vector<complex<float> > work(lwork);
3190  Vector<float> rwork(3*n);
3191  Conjugate(A); Conjugate(B);
3192  w.Reallocate(n);
3193  z.Reallocate(n,n);
3194  chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3195  z.GetData(), &n, work.GetData(), rwork.GetData(),
3196  &info.GetInfoRef());
3197 
3198 #ifdef SELDON_LAPACK_CHECK_INFO
3199  if (info.GetInfo() != 0)
3200  throw LapackError(info.GetInfo(), "GetEigenvalues",
3201  "Failed to find eigenvalues ");
3202 #endif
3203 
3204  Transpose(z);
3205  }
3206 
3207 
3208  template<class Prop1, class Prop2, class Allocator1,
3209  class Allocator2, class Allocator3>
3210  void GetEigenvalues(Matrix<complex<double>,
3211  Prop1, RowHermPacked, Allocator1>& A,
3212  Matrix<complex<double>,
3213  Prop2, RowHermPacked, Allocator2>& B,
3214  Vector<double, VectFull, Allocator3>& w,
3215  LapackInfo& info)
3216  {
3217 #ifdef SELDON_CHECK_DIMENSIONS
3218  if (A.GetM() != B.GetM())
3219  throw WrongDim("GetEigenvalues",
3220  "Matrix A and B must have the same size");
3221 #endif
3222 
3223  int itype = 1;
3224  int n = A.GetM();
3225  char uplo('L'), job('N');
3226  int lwork = 2*n;
3227  Vector<complex<double> > work(lwork);
3228  Vector<double> rwork(3*n);
3229  w.Reallocate(n);
3230  Conjugate(A); Conjugate(B);
3231  zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3232  A.GetData(), &n, work.GetData(), rwork.GetData(),
3233  &info.GetInfoRef());
3234 
3235 #ifdef SELDON_LAPACK_CHECK_INFO
3236  if (info.GetInfo() != 0)
3237  throw LapackError(info.GetInfo(), "GetEigenvalues",
3238  "Failed to find eigenvalues ");
3239 #endif
3240 
3241  }
3242 
3243 
3244  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3245  class Allocator3, class Allocator4>
3246  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
3247  Prop1, RowHermPacked, Allocator1>& A,
3248  Matrix<complex<double>,
3249  Prop2, RowHermPacked, Allocator2>& B,
3250  Vector<double, VectFull, Allocator3>& w,
3251  Matrix<complex<double>,
3252  General, RowMajor, Allocator4>& z,
3253  LapackInfo& info)
3254  {
3255 #ifdef SELDON_CHECK_DIMENSIONS
3256  if (A.GetM() != B.GetM())
3257  throw WrongDim("GetEigenvaluesEigenvectors",
3258  "Matrix A and B must have the same size");
3259 #endif
3260 
3261  int itype = 1;
3262  int n = A.GetM();
3263  char uplo('L'); char job('V'); int lwork = 2*n;
3264  Vector<complex<double> > work(lwork);
3265  Vector<double> rwork(3*n);
3266  w.Reallocate(n);
3267  z.Reallocate(n,n);
3268  Conjugate(A); Conjugate(B);
3269  zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3270  z.GetData(), &n, work.GetData(), rwork.GetData(),
3271  &info.GetInfoRef());
3272 
3273 #ifdef SELDON_LAPACK_CHECK_INFO
3274  if (info.GetInfo() != 0)
3275  throw LapackError(info.GetInfo(), "GetEigenvalues",
3276  "Failed to find eigenvalues ");
3277 #endif
3278 
3279  Transpose(z);
3280  }
3281 
3282 
3283  /* ColHermPacked */
3284 
3285 
3286  template<class Prop1, class Prop2, class Allocator1,
3287  class Allocator2, class Allocator3>
3288  void GetEigenvalues(Matrix<complex<float>,
3289  Prop1, ColHermPacked, Allocator1>& A,
3290  Matrix<complex<float>,
3291  Prop2, ColHermPacked, Allocator2>& B,
3292  Vector<float, VectFull, Allocator3>& w,
3293  LapackInfo& info)
3294  {
3295 #ifdef SELDON_CHECK_DIMENSIONS
3296  if (A.GetM() != B.GetM())
3297  throw WrongDim("GetEigenvalues",
3298  "Matrix A and B must have the same size");
3299 #endif
3300 
3301  int itype = 1;
3302  int n = A.GetM();
3303  char uplo('U'); char job('N');
3304  int lwork = 2*n;
3305  Vector<complex<float> > work(lwork);
3306  Vector<float> rwork(3*n);
3307  w.Reallocate(n);
3308  chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(),
3309  w.GetData(), A.GetData(), &n, work.GetData(),
3310  rwork.GetData(), &info.GetInfoRef());
3311 
3312 #ifdef SELDON_LAPACK_CHECK_INFO
3313  if (info.GetInfo() != 0)
3314  throw LapackError(info.GetInfo(), "GetEigenvalues",
3315  "Failed to find eigenvalues ");
3316 #endif
3317 
3318  }
3319 
3320 
3321  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3322  class Allocator3, class Allocator4>
3323  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
3324  Prop1, ColHermPacked, Allocator1>& A,
3325  Matrix<complex<float>,
3326  Prop2, ColHermPacked, Allocator2>& B,
3327  Vector<float, VectFull, Allocator3>& w,
3328  Matrix<complex<float>,
3329  General, ColMajor, Allocator4>& z,
3330  LapackInfo& info)
3331  {
3332 #ifdef SELDON_CHECK_DIMENSIONS
3333  if (A.GetM() != B.GetM())
3334  throw WrongDim("GetEigenvaluesEigenvectors",
3335  "Matrix A and B must have the same size");
3336 #endif
3337 
3338  int itype = 1;
3339  int n = A.GetM();
3340  char uplo('U'); char job('V');
3341  int lwork = 3*n;
3342  Vector<complex<float> > work(lwork);
3343  Vector<float> rwork(3*n);
3344  w.Reallocate(n);
3345  z.Reallocate(n,n);
3346  chpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3347  z.GetData(), &n, work.GetData(), rwork.GetData(),
3348  &info.GetInfoRef());
3349 
3350 #ifdef SELDON_LAPACK_CHECK_INFO
3351  if (info.GetInfo() != 0)
3352  throw LapackError(info.GetInfo(), "GetEigenvalues",
3353  "Failed to find eigenvalues ");
3354 #endif
3355 
3356  }
3357 
3358 
3359  template<class Prop1, class Prop2, class Allocator1,
3360  class Allocator2, class Allocator3>
3361  void GetEigenvalues(Matrix<complex<double>,
3362  Prop1, ColHermPacked, Allocator1>& A,
3363  Matrix<complex<double>,
3364  Prop2, ColHermPacked, Allocator2>& B,
3365  Vector<double, VectFull, Allocator3>& w,
3366  LapackInfo& info)
3367  {
3368 #ifdef SELDON_CHECK_DIMENSIONS
3369  if (A.GetM() != B.GetM())
3370  throw WrongDim("GetEigenvalues",
3371  "Matrix A and B must have the same size");
3372 #endif
3373 
3374  int itype = 1;
3375  int n = A.GetM();
3376  char uplo('U'); char job('N');
3377  int lwork = 2*n;
3378  Vector<complex<double> > work(lwork);
3379  Vector<double> rwork(3*n);
3380  w.Reallocate(n);
3381  zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(),
3382  w.GetData(), A.GetData(), &n, work.GetData(),
3383  rwork.GetData(), &info.GetInfoRef());
3384 
3385 #ifdef SELDON_LAPACK_CHECK_INFO
3386  if (info.GetInfo() != 0)
3387  throw LapackError(info.GetInfo(), "GetEigenvalues",
3388  "Failed to find eigenvalues ");
3389 #endif
3390 
3391  }
3392 
3393 
3394  template<class Prop1, class Prop2, class Allocator1, class Allocator2,
3395  class Allocator3, class Allocator4>
3396  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
3397  Prop1, ColHermPacked, Allocator1>& A,
3398  Matrix<complex<double>,
3399  Prop2, ColHermPacked, Allocator2>& B,
3400  Vector<double, VectFull, Allocator3>& w,
3401  Matrix<complex<double>,
3402  General, ColMajor, Allocator4>& z,
3403  LapackInfo& info)
3404  {
3405 #ifdef SELDON_CHECK_DIMENSIONS
3406  if (A.GetM() != B.GetM())
3407  throw WrongDim("GetEigenvaluesEigenvectors",
3408  "Matrix A and B must have the same size");
3409 #endif
3410 
3411  int itype = 1;
3412  int n = A.GetM();
3413  char uplo('U'); char job('V');
3414  int lwork = 3*n;
3415  Vector<complex<double> > work(lwork);
3416  Vector<double> rwork(3*n);
3417  w.Reallocate(n);
3418  z.Reallocate(n,n);
3419  zhpgv_(&itype, &job, &uplo, &n, A.GetData(), B.GetData(), w.GetData(),
3420  z.GetData(), &n, work.GetData(), rwork.GetData(),
3421  &info.GetInfoRef());
3422 
3423 #ifdef SELDON_LAPACK_CHECK_INFO
3424  if (info.GetInfo() != 0)
3425  throw LapackError(info.GetInfo(), "GetEigenvalues",
3426  "Failed to find eigenvalues ");
3427 #endif
3428 
3429  }
3430 
3431 
3432  /* RowMajor */
3433 
3434 
3435  template<class Prop1, class Prop2, class Allocator1,
3436  class Allocator2, class Allocator3,
3437  class Allocator4, class Allocator5>
3438  void GetEigenvalues(Matrix<float, Prop1, RowMajor, Allocator1>& A,
3439  Matrix<float, Prop2, RowMajor, Allocator2>& B,
3440  Vector<float, VectFull, Allocator3>& alpha_real,
3441  Vector<float, VectFull, Allocator4>& alpha_imag,
3442  Vector<float, VectFull, Allocator5>& beta,
3443  LapackInfo& info)
3444  {
3445 #ifdef SELDON_CHECK_DIMENSIONS
3446  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3447  throw WrongDim("GetEigenvalues",
3448  "Matrix A and B must be squared");
3449 
3450  if (A.GetM() != B.GetM())
3451  throw WrongDim("GetEigenvalues",
3452  "Matrix A and B must have the same size");
3453 #endif
3454 
3455  int n = A.GetM();
3456  char jobvl('N'), jobvr('N');
3457  int lwork = 8*n+16; Vector<float> work(lwork);
3458  alpha_real.Reallocate(n);
3459  alpha_imag.Reallocate(n);
3460  beta.Reallocate(n);
3461  sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3462  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
3463  A.GetData(), &n, A.GetData(), &n,
3464  work.GetData(), &lwork, &info.GetInfoRef());
3465 
3466 #ifdef SELDON_LAPACK_CHECK_INFO
3467  if (info.GetInfo() != 0)
3468  throw LapackError(info.GetInfo(), "GetEigenvalues",
3469  "Failed to find eigenvalues ");
3470 #endif
3471 
3472  }
3473 
3474 
3475  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3476  class Allocator2, class Allocator3, class Allocator4,
3477  class Allocator5, class Allocator6>
3478  void GetEigenvaluesEigenvectors(Matrix<float, Prop1, RowMajor, Allocator1>& A,
3479  Matrix<float, Prop2, RowMajor, Allocator2>& B,
3480  Vector<float, VectFull, Allocator3>& alpha_real,
3481  Vector<float, VectFull, Allocator4>& alpha_imag,
3482  Vector<float, VectFull, Allocator5>& beta,
3483  Matrix<float, Prop3, RowMajor, Allocator6>& V,
3484  LapackInfo& info)
3485  {
3486 #ifdef SELDON_CHECK_DIMENSIONS
3487  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3488  throw WrongDim("GetEigenvaluesEigenvectors",
3489  "Matrix A and B must be squared");
3490 
3491  if (A.GetM() != B.GetM())
3492  throw WrongDim("GetEigenvaluesEigenvectors",
3493  "Matrix A and B must have the same size");
3494 #endif
3495 
3496  int n = A.GetM();
3497  char jobvl('V'), jobvr('N');
3498  int lwork = 8*n+16; Vector<float> work(lwork);
3499  alpha_real.Reallocate(n);
3500  alpha_imag.Reallocate(n);
3501  beta.Reallocate(n);
3502  V.Reallocate(n,n);
3503  sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3504  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
3505  V.GetData(), &n, V.GetData(), &n,
3506  work.GetData(), &lwork, &info.GetInfoRef());
3507 
3508 
3509 #ifdef SELDON_LAPACK_CHECK_INFO
3510  if (info.GetInfo() != 0)
3511  throw LapackError(info.GetInfo(), "GetEigenvalues",
3512  "Failed to find eigenvalues ");
3513 #endif
3514 
3515  Transpose(V);
3516  // conjugate if necessary
3517  int i = 0;
3518  while (i < n)
3519  {
3520  if (i < (n-1))
3521  if (alpha_real(i) == alpha_real(i+1))
3522  {
3523  for (int j = 0; j < n; j++)
3524  V(j,i+1) = -V(j,i+1);
3525 
3526  i++;
3527  }
3528 
3529  i++;
3530  }
3531  }
3532 
3533 
3534  template<class Prop1, class Prop2, class Allocator1,
3535  class Allocator2, class Allocator4, class Allocator5>
3536  void GetEigenvalues(Matrix<complex<float>, Prop1, RowMajor, Allocator1>& A,
3537  Matrix<complex<float>, Prop2, RowMajor, Allocator2>& B,
3538  Vector<complex<float>, VectFull, Allocator4>& alpha,
3539  Vector<complex<float>, VectFull, Allocator5>& beta,
3540  LapackInfo& info)
3541  {
3542 #ifdef SELDON_CHECK_DIMENSIONS
3543  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3544  throw WrongDim("GetEigenvalues",
3545  "Matrix A and B must be squared");
3546 
3547  if (A.GetM() != B.GetM())
3548  throw WrongDim("GetEigenvalues",
3549  "Matrix A and B must have the same size");
3550 #endif
3551 
3552  int n = A.GetM();
3553  char jobvl('N'), jobvr('N'); int lwork = 2*n;
3554  Vector<complex<float> > work(lwork);
3555  Vector<float> rwork(8*n);
3556  alpha.Reallocate(n);
3557  beta.Reallocate(n);
3558  cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3559  alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n,
3560  work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef());
3561 
3562 #ifdef SELDON_LAPACK_CHECK_INFO
3563  if (info.GetInfo() != 0)
3564  throw LapackError(info.GetInfo(), "GetEigenvalues",
3565  "Failed to find eigenvalues ");
3566 #endif
3567 
3568  }
3569 
3570 
3571  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3572  class Allocator2, class Allocator4,
3573  class Allocator5, class Allocator6>
3574  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
3575  Prop1, RowMajor, Allocator1>& A,
3576  Matrix<complex<float>,
3577  Prop2, RowMajor, Allocator2>& B,
3578  Vector<complex<float>,
3579  VectFull, Allocator4>& alpha,
3580  Vector<complex<float>,
3581  VectFull, Allocator5>& beta,
3582  Matrix<complex<float>,
3583  Prop3, RowMajor, Allocator6>& V,
3584  LapackInfo& info)
3585  {
3586 #ifdef SELDON_CHECK_DIMENSIONS
3587  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3588  throw WrongDim("GetEigenvaluesEigenvectors",
3589  "Matrix A and B must be squared");
3590 
3591  if (A.GetM() != B.GetM())
3592  throw WrongDim("GetEigenvaluesEigenvectors",
3593  "Matrix A and B must have the same size");
3594 #endif
3595 
3596  int n = A.GetM();
3597  char jobvl('V'), jobvr('N');
3598  int lwork = 2*n; Vector<complex<float> > work(lwork);
3599  Vector<float> rwork(8*n);
3600  alpha.Reallocate(n);
3601  beta.Reallocate(n);
3602  V.Reallocate(n,n);
3603  cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, alpha.GetData(),
3604  beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(),
3605  &lwork, rwork.GetData(), &info.GetInfoRef());
3606 
3607 
3608 #ifdef SELDON_LAPACK_CHECK_INFO
3609  if (info.GetInfo() != 0)
3610  throw LapackError(info.GetInfo(), "GetEigenvalues",
3611  "Failed to find eigenvalues ");
3612 #endif
3613 
3614  TransposeConj(V);
3615  }
3616 
3617 
3618  template<class Prop1, class Prop2, class Allocator1,
3619  class Allocator2, class Allocator3,
3620  class Allocator4, class Allocator5>
3621  void GetEigenvalues(Matrix<double, Prop1, RowMajor, Allocator1>& A,
3622  Matrix<double, Prop2, RowMajor, Allocator2>& B,
3623  Vector<double, VectFull, Allocator3>& alpha_real,
3624  Vector<double, VectFull, Allocator4>& alpha_imag,
3625  Vector<double, VectFull, Allocator5>& beta,
3626  LapackInfo& info)
3627  {
3628 #ifdef SELDON_CHECK_DIMENSIONS
3629  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3630  throw WrongDim("GetEigenvalues",
3631  "Matrix A and B must be squared");
3632 
3633  if (A.GetM() != B.GetM())
3634  throw WrongDim("GetEigenvalues",
3635  "Matrix A and B must have the same size");
3636 #endif
3637 
3638  int n = A.GetM();
3639  char jobvl('N'), jobvr('N');
3640  int lwork = 8*n+16; Vector<double> work(lwork);
3641  alpha_real.Reallocate(n);
3642  alpha_imag.Reallocate(n);
3643  beta.Reallocate(n);
3644  dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3645  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
3646  A.GetData(), &n, A.GetData(), &n,
3647  work.GetData(), &lwork, &info.GetInfoRef());
3648 
3649 #ifdef SELDON_LAPACK_CHECK_INFO
3650  if (info.GetInfo() != 0)
3651  throw LapackError(info.GetInfo(), "GetEigenvalues",
3652  "Failed to find eigenvalues ");
3653 #endif
3654 
3655  }
3656 
3657 
3658  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3659  class Allocator2, class Allocator3, class Allocator4,
3660  class Allocator5, class Allocator6>
3661  void GetEigenvaluesEigenvectors(Matrix<double, Prop1, RowMajor, Allocator1>& A,
3662  Matrix<double, Prop2, RowMajor, Allocator2>& B,
3663  Vector<double, VectFull, Allocator3>& alpha_real,
3664  Vector<double, VectFull, Allocator4>& alpha_imag,
3665  Vector<double, VectFull, Allocator5>& beta,
3666  Matrix<double, Prop3, RowMajor, Allocator6>& V,
3667  LapackInfo& info)
3668  {
3669 #ifdef SELDON_CHECK_DIMENSIONS
3670  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3671  throw WrongDim("GetEigenvaluesEigenvectors",
3672  "Matrix A and B must be squared");
3673 
3674  if (A.GetM() != B.GetM())
3675  throw WrongDim("GetEigenvaluesEigenvectors",
3676  "Matrix A and B must have the same size");
3677 #endif
3678 
3679  int n = A.GetM();
3680  char jobvl('V'), jobvr('N');
3681  int lwork = 8*n+16; Vector<double> work(lwork);
3682  alpha_real.Reallocate(n);
3683  alpha_imag.Reallocate(n);
3684  beta.Reallocate(n);
3685  V.Reallocate(n, n);
3686  dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3687  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
3688  V.GetData(), &n, V.GetData(), &n,
3689  work.GetData(), &lwork, &info.GetInfoRef());
3690 
3691 #ifdef SELDON_LAPACK_CHECK_INFO
3692  if (info.GetInfo() != 0)
3693  throw LapackError(info.GetInfo(), "GetEigenvalues",
3694  "Failed to find eigenvalues ");
3695 #endif
3696 
3697  Transpose(V);
3698  // conjugate if necessary
3699  int i = 0;
3700  while (i < n)
3701  {
3702  if (i < (n-1))
3703  if (alpha_real(i) == alpha_real(i+1))
3704  {
3705  for (int j = 0; j < n; j++)
3706  V(j,i+1) = -V(j,i+1);
3707 
3708  i++;
3709  }
3710 
3711  i++;
3712  }
3713  }
3714 
3715 
3716  template<class Prop1, class Prop2, class Allocator1,
3717  class Allocator2, class Allocator4, class Allocator5>
3718  void GetEigenvalues(Matrix<complex<double>, Prop1, RowMajor, Allocator1>& A,
3719  Matrix<complex<double>, Prop2, RowMajor, Allocator2>& B,
3720  Vector<complex<double>, VectFull, Allocator4>& alpha,
3721  Vector<complex<double>, VectFull, Allocator5>& beta,
3722  LapackInfo& info)
3723  {
3724 #ifdef SELDON_CHECK_DIMENSIONS
3725  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3726  throw WrongDim("GetEigenvalues",
3727  "Matrix A and B must be squared");
3728 
3729  if (A.GetM() != B.GetM())
3730  throw WrongDim("GetEigenvalues",
3731  "Matrix A and B must have the same size");
3732 #endif
3733 
3734  int n = A.GetM();
3735  char jobvl('N'), jobvr('N'); int lwork = 2*n;
3736  Vector<complex<double> > work(lwork);
3737  Vector<double> rwork(8*n);
3738  alpha.Reallocate(n);
3739  beta.Reallocate(n);
3740  zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3741  alpha.GetData(), beta.GetData(), A.GetData(), &n, A.GetData(), &n,
3742  work.GetData(), &lwork, rwork.GetData(), &info.GetInfoRef());
3743 
3744 #ifdef SELDON_LAPACK_CHECK_INFO
3745  if (info.GetInfo() != 0)
3746  throw LapackError(info.GetInfo(), "GetEigenvalues",
3747  "Failed to find eigenvalues ");
3748 #endif
3749 
3750  }
3751 
3752 
3753  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3754  class Allocator2, class Allocator4,
3755  class Allocator5, class Allocator6>
3756  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
3757  Prop1, RowMajor, Allocator1>& A,
3758  Matrix<complex<double>,
3759  Prop2, RowMajor, Allocator2>& B,
3760  Vector<complex<double>,
3761  VectFull, Allocator4>& alpha,
3762  Vector<complex<double>,
3763  VectFull, Allocator5>& beta,
3764  Matrix<complex<double>,
3765  Prop3, RowMajor, Allocator6>& V,
3766  LapackInfo& info)
3767  {
3768 #ifdef SELDON_CHECK_DIMENSIONS
3769  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3770  throw WrongDim("GetEigenvaluesEigenvectors",
3771  "Matrix A and B must be squared");
3772 
3773  if (A.GetM() != B.GetM())
3774  throw WrongDim("GetEigenvaluesEigenvectors",
3775  "Matrix A and B must have the same size");
3776 #endif
3777 
3778  int n = A.GetM();
3779  char jobvl('V'), jobvr('N');
3780  int lwork = 2*n; Vector<complex<double> > work(lwork);
3781  Vector<double> rwork(8*n);
3782  alpha.Reallocate(n);
3783  beta.Reallocate(n);
3784  V.Reallocate(n,n);
3785  zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n, alpha.GetData(),
3786  beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(),
3787  &lwork, rwork.GetData(), &info.GetInfoRef());
3788 
3789 #ifdef SELDON_LAPACK_CHECK_INFO
3790  if (info.GetInfo() != 0)
3791  throw LapackError(info.GetInfo(), "GetEigenvalues",
3792  "Failed to find eigenvalues ");
3793 #endif
3794 
3795  TransposeConj(V);
3796  }
3797 
3798 
3799  /* ColMajor */
3800 
3801 
3802  template<class Prop1, class Prop2, class Allocator1,
3803  class Allocator2, class Allocator3,
3804  class Allocator4, class Allocator5>
3805  void GetEigenvalues(Matrix<float, Prop1, ColMajor, Allocator1>& A,
3806  Matrix<float, Prop2, ColMajor, Allocator2>& B,
3807  Vector<float, VectFull, Allocator3>& alpha_real,
3808  Vector<float, VectFull, Allocator4>& alpha_imag,
3809  Vector<float, VectFull, Allocator5>& beta,
3810  LapackInfo& info)
3811  {
3812 #ifdef SELDON_CHECK_DIMENSIONS
3813  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3814  throw WrongDim("GetEigenvalues",
3815  "Matrix A and B must be squared");
3816 
3817  if (A.GetM() != B.GetM())
3818  throw WrongDim("GetEigenvalues",
3819  "Matrix A and B must have the same size");
3820 #endif
3821 
3822  int n = A.GetM();
3823  char jobvl('N'), jobvr('N');
3824  int lwork = 8*n+16; Vector<float> work(lwork);
3825  alpha_real.Reallocate(n);
3826  alpha_imag.Reallocate(n);
3827  beta.Reallocate(n);
3828 
3829  sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
3830  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
3831  A.GetData(), &n, A.GetData(), &n,
3832  work.GetData(), &lwork, &info.GetInfoRef());
3833 
3834 #ifdef SELDON_LAPACK_CHECK_INFO
3835  if (info.GetInfo() != 0)
3836  throw LapackError(info.GetInfo(), "GetEigenvalues",
3837  "Failed to find eigenvalues ");
3838 #endif
3839 
3840  }
3841 
3842 
3843  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3844  class Allocator2, class Allocator3, class Allocator4,
3845  class Allocator5, class Allocator6>
3846  void GetEigenvaluesEigenvectors(Matrix<float, Prop1, ColMajor, Allocator1>& A,
3847  Matrix<float, Prop2, ColMajor, Allocator2>& B,
3848  Vector<float,
3849  VectFull, Allocator3>& alpha_real,
3850  Vector<float,
3851  VectFull, Allocator4>& alpha_imag,
3852  Vector<float, VectFull, Allocator5>& beta,
3853  Matrix<float, Prop3, ColMajor, Allocator6>& V,
3854  LapackInfo& info)
3855  {
3856 #ifdef SELDON_CHECK_DIMENSIONS
3857  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3858  throw WrongDim("GetEigenvaluesEigenvectors",
3859  "Matrix A and B must be squared");
3860 
3861  if (A.GetM() != B.GetM())
3862  throw WrongDim("GetEigenvaluesEigenvectors",
3863  "Matrix A and B must have the same size");
3864 #endif
3865 
3866  int n = A.GetM();
3867  char jobvl('V'), jobvr('N');
3868  int lwork = 8*n+16; Vector<float> work(lwork);
3869  alpha_real.Reallocate(n);
3870  alpha_imag.Reallocate(n);
3871  beta.Reallocate(n);
3872  V.Reallocate(n,n);
3873  sggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
3874  &n, alpha_real.GetData(), alpha_imag.GetData(),
3875  beta.GetData(), V.GetData(), &n, V.GetData(), &n,
3876  work.GetData(), &lwork, &info.GetInfoRef());
3877 
3878 #ifdef SELDON_LAPACK_CHECK_INFO
3879  if (info.GetInfo() != 0)
3880  throw LapackError(info.GetInfo(), "GetEigenvalues",
3881  "Failed to find eigenvalues ");
3882 #endif
3883 
3884  }
3885 
3886 
3887  template<class Prop1, class Prop2, class Allocator1,
3888  class Allocator2, class Allocator4, class Allocator5>
3889  void GetEigenvalues(Matrix<complex<float>, Prop1, ColMajor, Allocator1>& A,
3890  Matrix<complex<float>, Prop2, ColMajor, Allocator2>& B,
3891  Vector<complex<float>, VectFull, Allocator4>& alpha,
3892  Vector<complex<float>, VectFull, Allocator5>& beta,
3893  LapackInfo& info)
3894  {
3895 #ifdef SELDON_CHECK_DIMENSIONS
3896  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3897  throw WrongDim("GetEigenvalues",
3898  "Matrix A and B must be squared");
3899 
3900  if (A.GetM() != B.GetM())
3901  throw WrongDim("GetEigenvalues",
3902  "Matrix A and B must have the same size");
3903 #endif
3904 
3905  int n = A.GetM();
3906  char jobvl('N'), jobvr('N'); int lwork = 2*n;
3907  Vector<complex<float> > work(lwork);
3908  Vector<float> rwork(8*n);
3909  alpha.Reallocate(n);
3910  beta.Reallocate(n);
3911 
3912  cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
3913  &n, alpha.GetData(),
3914  beta.GetData(), A.GetData(), &n, A.GetData(), &n, work.GetData(),
3915  &lwork, rwork.GetData(), &info.GetInfoRef());
3916 
3917 #ifdef SELDON_LAPACK_CHECK_INFO
3918  if (info.GetInfo() != 0)
3919  throw LapackError(info.GetInfo(), "GetEigenvalues",
3920  "Failed to find eigenvalues ");
3921 #endif
3922 
3923  }
3924 
3925 
3926  template<class Prop1, class Prop2, class Prop3, class Allocator1,
3927  class Allocator2, class Allocator4,
3928  class Allocator5, class Allocator6>
3929  void GetEigenvaluesEigenvectors(Matrix<complex<float>,
3930  Prop1, ColMajor, Allocator1>& A,
3931  Matrix<complex<float>,
3932  Prop2, ColMajor, Allocator2>& B,
3933  Vector<complex<float>,
3934  VectFull, Allocator4>& alpha,
3935  Vector<complex<float>,
3936  VectFull, Allocator5>& beta,
3937  Matrix<complex<float>,
3938  Prop3, ColMajor, Allocator6>& V,
3939  LapackInfo& info)
3940  {
3941 #ifdef SELDON_CHECK_DIMENSIONS
3942  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3943  throw WrongDim("GetEigenvaluesEigenvectors",
3944  "Matrix A and B must be squared");
3945 
3946  if (A.GetM() != B.GetM())
3947  throw WrongDim("GetEigenvaluesEigenvectors",
3948  "Matrix A and B must have the same size");
3949 #endif
3950 
3951  int n = A.GetM();
3952  char jobvl('N'), jobvr('V'); int lwork = 2*n;
3953  Vector<complex<float> > work(lwork);
3954  Vector<float> rwork(8*n);
3955  alpha.Reallocate(n);
3956  beta.Reallocate(n);
3957  V.Reallocate(n,n);
3958  cggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
3959  &n, alpha.GetData(),
3960  beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(),
3961  &lwork, rwork.GetData(), &info.GetInfoRef());
3962 
3963 #ifdef SELDON_LAPACK_CHECK_INFO
3964  if (info.GetInfo() != 0)
3965  throw LapackError(info.GetInfo(), "GetEigenvalues",
3966  "Failed to find eigenvalues ");
3967 #endif
3968 
3969  }
3970 
3971 
3972  template<class Prop1, class Prop2, class Allocator1,
3973  class Allocator2, class Allocator3,
3974  class Allocator4, class Allocator5>
3975  void GetEigenvalues(Matrix<double, Prop1, ColMajor, Allocator1>& A,
3976  Matrix<double, Prop2, ColMajor, Allocator2>& B,
3977  Vector<double, VectFull, Allocator3>& alpha_real,
3978  Vector<double, VectFull, Allocator4>& alpha_imag,
3979  Vector<double, VectFull, Allocator5>& beta,
3980  LapackInfo& info)
3981  {
3982 #ifdef SELDON_CHECK_DIMENSIONS
3983  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
3984  throw WrongDim("GetEigenvalues",
3985  "Matrix A and B must be squared");
3986 
3987  if (A.GetM() != B.GetM())
3988  throw WrongDim("GetEigenvalues",
3989  "Matrix A and B must have the same size");
3990 #endif
3991 
3992  int n = A.GetM();
3993  char jobvl('N'), jobvr('N');
3994  int lwork = 8*n+16; Vector<double> work(lwork);
3995  alpha_real.Reallocate(n);
3996  alpha_imag.Reallocate(n);
3997  beta.Reallocate(n);
3998 
3999  dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(), &n,
4000  alpha_real.GetData(), alpha_imag.GetData(), beta.GetData(),
4001  A.GetData(), &n, A.GetData(), &n,
4002  work.GetData(), &lwork, &info.GetInfoRef());
4003 
4004 #ifdef SELDON_LAPACK_CHECK_INFO
4005  if (info.GetInfo() != 0)
4006  throw LapackError(info.GetInfo(), "GetEigenvalues",
4007  "Failed to find eigenvalues ");
4008 #endif
4009 
4010  }
4011 
4012 
4013  template<class Prop1, class Prop2, class Prop3, class Allocator1,
4014  class Allocator2, class Allocator3, class Allocator4,
4015  class Allocator5, class Allocator6>
4016  void GetEigenvaluesEigenvectors(Matrix<double, Prop1, ColMajor, Allocator1>& A,
4017  Matrix<double, Prop2, ColMajor, Allocator2>& B,
4018  Vector<double,
4019  VectFull, Allocator3>& alpha_real,
4020  Vector<double,
4021  VectFull, Allocator4>& alpha_imag,
4022  Vector<double, VectFull, Allocator5>& beta,
4023  Matrix<double, Prop3, ColMajor, Allocator6>& V,
4024  LapackInfo& info)
4025  {
4026 #ifdef SELDON_CHECK_DIMENSIONS
4027  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
4028  throw WrongDim("GetEigenvaluesEigenvectors",
4029  "Matrix A and B must be squared");
4030 
4031  if (A.GetM() != B.GetM())
4032  throw WrongDim("GetEigenvaluesEigenvectors",
4033  "Matrix A and B must have the same size");
4034 #endif
4035 
4036  int n = A.GetM();
4037  char jobvl('V'), jobvr('N');
4038  int lwork = 8*n+16; Vector<double> work(lwork);
4039  alpha_real.Reallocate(n);
4040  alpha_imag.Reallocate(n);
4041  beta.Reallocate(n);
4042  V.Reallocate(n,n);
4043 
4044  dggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
4045  &n, alpha_real.GetData(), alpha_imag.GetData(),
4046  beta.GetData(), V.GetData(), &n, V.GetData(), &n,
4047  work.GetData(), &lwork, &info.GetInfoRef());
4048 
4049 #ifdef SELDON_LAPACK_CHECK_INFO
4050  if (info.GetInfo() != 0)
4051  throw LapackError(info.GetInfo(), "GetEigenvalues",
4052  "Failed to find eigenvalues ");
4053 #endif
4054 
4055  }
4056 
4057 
4058  template<class Prop1, class Prop2, class Allocator1,
4059  class Allocator2, class Allocator4, class Allocator5>
4060  void GetEigenvalues(Matrix<complex<double>, Prop1, ColMajor, Allocator1>& A,
4061  Matrix<complex<double>, Prop2, ColMajor, Allocator2>& B,
4062  Vector<complex<double>, VectFull, Allocator4>& alpha,
4063  Vector<complex<double>, VectFull, Allocator5>& beta,
4064  LapackInfo& info)
4065  {
4066 #ifdef SELDON_CHECK_DIMENSIONS
4067  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
4068  throw WrongDim("GetEigenvalues",
4069  "Matrix A and B must be squared");
4070 
4071  if (A.GetM() != B.GetM())
4072  throw WrongDim("GetEigenvalues",
4073  "Matrix A and B must have the same size");
4074 #endif
4075 
4076  int n = A.GetM();
4077  char jobvl('N'), jobvr('N'); int lwork = 2*n;
4078  Vector<complex<double> > work(lwork);
4079  Vector<double> rwork(8*n);
4080  alpha.Reallocate(n);
4081  beta.Reallocate(n);
4082 
4083  zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
4084  &n, alpha.GetData(),
4085  beta.GetData(), A.GetData(), &n, A.GetData(), &n, work.GetData(),
4086  &lwork, rwork.GetData(), &info.GetInfoRef());
4087 
4088 #ifdef SELDON_LAPACK_CHECK_INFO
4089  if (info.GetInfo() != 0)
4090  throw LapackError(info.GetInfo(), "GetEigenvalues",
4091  "Failed to find eigenvalues ");
4092 #endif
4093 
4094  }
4095 
4096 
4097  template<class Prop1, class Prop2, class Prop3, class Allocator1,
4098  class Allocator2, class Allocator4,
4099  class Allocator5, class Allocator6>
4100  void GetEigenvaluesEigenvectors(Matrix<complex<double>,
4101  Prop1, ColMajor, Allocator1>& A,
4102  Matrix<complex<double>,
4103  Prop2, ColMajor, Allocator2>& B,
4104  Vector<complex<double>,
4105  VectFull, Allocator4>& alpha,
4106  Vector<complex<double>,
4107  VectFull, Allocator5>& beta,
4108  Matrix<complex<double>,
4109  Prop3, ColMajor, Allocator6>& V,
4110  LapackInfo& info)
4111  {
4112 #ifdef SELDON_CHECK_DIMENSIONS
4113  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()) )
4114  throw WrongDim("GetEigenvaluesEigenvectors",
4115  "Matrix A and B must be squared");
4116 
4117  if (A.GetM() != B.GetM())
4118  throw WrongDim("GetEigenvaluesEigenvectors",
4119  "Matrix A and B must have the same size");
4120 #endif
4121 
4122  int n = A.GetM();
4123  char jobvl('N'), jobvr('V'); int lwork = 2*n;
4124  Vector<complex<double> > work(lwork);
4125  Vector<double> rwork(8*n);
4126  alpha.Reallocate(n);
4127  beta.Reallocate(n);
4128  V.Reallocate(n,n);
4129  zggev_(&jobvl, &jobvr, &n, A.GetData(), &n, B.GetData(),
4130  &n, alpha.GetData(),
4131  beta.GetData(), V.GetData(), &n, V.GetData(), &n, work.GetData(),
4132  &lwork, rwork.GetData(), &info.GetInfoRef());
4133 
4134 #ifdef SELDON_LAPACK_CHECK_INFO
4135  if (info.GetInfo() != 0)
4136  throw LapackError(info.GetInfo(), "GetEigenvalues",
4137  "Failed to find eigenvalues ");
4138 #endif
4139 
4140  }
4141 
4142 
4143  // GENERALIZED EIGENVALUE PROBLEM //
4145 
4146 
4148  // SINGULAR VALUE DECOMPOSITION //
4149 
4150 
4151  /* RowMajor */
4152 
4153 
4154  template<class Prop1, class Allocator1, class Allocator2,
4155  class Allocator3, class Allocator4>
4156  void GetSVD(Matrix<float, Prop1, RowMajor, Allocator1>& A,
4157  Vector<float, VectFull, Allocator4>& lambda,
4158  Matrix<float, General, RowMajor, Allocator2>& u,
4159  Matrix<float, General, RowMajor, Allocator3>& v,
4160  LapackInfo& info)
4161  {
4162  int m = A.GetM();
4163  int n = A.GetN();
4164  char jobl('A'), jobr('A');
4165  int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n));
4166  Vector<float> work(lwork);
4167  lambda.Reallocate(min(m, n));
4168  lambda.Zero();
4169  u.Reallocate(m, m);
4170  u.Zero();
4171  v.Reallocate(n, n);
4172  v.Zero();
4173  sgesvd_(&jobl, &jobr, &n, &m, A.GetData(), &n, lambda.GetData(),
4174  v.GetData(), &n, u.GetData(), &m, work.GetData(),
4175  &lwork, &info.GetInfoRef());
4176 
4177 #ifdef SELDON_LAPACK_CHECK_INFO
4178  if (info.GetInfo() != 0)
4179  throw LapackError(info.GetInfo(), "GetSVD",
4180  "Failed to find singular value decomposition");
4181 #endif
4182 
4183  }
4184 
4185 
4186  template<class Prop1, class Allocator1, class Allocator2,
4187  class Allocator3, class Allocator4>
4188  void GetSVD(Matrix<complex<float>, Prop1, RowMajor, Allocator1>& A,
4189  Vector<float, VectFull, Allocator4>& lambda,
4190  Matrix<complex<float>, General, RowMajor, Allocator2>& u,
4191  Matrix<complex<float>, General, RowMajor, Allocator3>& v,
4192  LapackInfo& info)
4193  {
4194  int m = A.GetM();
4195  int n = A.GetN();
4196  char jobl('A'), jobr('A');
4197  int lwork = 2 * min(m, n) + max(m, n);
4198  Vector<complex<float> > work(lwork);
4199  Vector<float> rwork(5 * min(m, n));
4200  lambda.Reallocate(min(m, n));
4201  lambda.Zero();
4202  u.Reallocate(m, m);
4203  u.Zero();
4204  v.Reallocate(n, n);
4205  v.Zero();
4206  cgesvd_(&jobl, &jobr, &n, &m, A.GetDataVoid(), &n, lambda.GetData(),
4207  v.GetDataVoid(), &n, u.GetDataVoid(), &m, work.GetDataVoid(),
4208  &lwork, rwork.GetData(), &info.GetInfoRef());
4209 
4210 #ifdef SELDON_LAPACK_CHECK_INFO
4211  if (info.GetInfo() != 0)
4212  throw LapackError(info.GetInfo(), "GetSVD",
4213  "Failed to find singular value decomposition");
4214 #endif
4215 
4216  }
4217 
4218 
4219  template<class Prop1, class Allocator1, class Allocator2,
4220  class Allocator3, class Allocator4>
4221  void GetSVD(Matrix<double, Prop1, RowMajor, Allocator1>& A,
4222  Vector<double, VectFull, Allocator4>& lambda,
4223  Matrix<double, General, RowMajor, Allocator2>& u,
4224  Matrix<double, General, RowMajor, Allocator3>& v,
4225  LapackInfo& info)
4226  {
4227  int m = A.GetM();
4228  int n = A.GetN();
4229  char jobl('A'), jobr('A');
4230  int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n));
4231  Vector<double> work(lwork);
4232  lambda.Reallocate(min(m, n));
4233  lambda.Zero();
4234  u.Reallocate(m, m);
4235  u.Zero();
4236  v.Reallocate(n, n);
4237  v.Zero();
4238  dgesvd_(&jobl, &jobr, &n, &m, A.GetData(), &n, lambda.GetData(),
4239  v.GetData(), &n, u.GetData(), &m, work.GetData(),
4240  &lwork, &info.GetInfoRef());
4241 
4242 #ifdef SELDON_LAPACK_CHECK_INFO
4243  if (info.GetInfo() != 0)
4244  throw LapackError(info.GetInfo(), "GetSVD",
4245  "Failed to find singular value decomposition");
4246 #endif
4247 
4248  }
4249 
4250 
4251  template<class Prop1, class Allocator1, class Allocator2,
4252  class Allocator3, class Allocator4>
4253  void GetSVD(Matrix<complex<double>, Prop1, RowMajor, Allocator1>& A,
4254  Vector<double, VectFull, Allocator4>& lambda,
4255  Matrix<complex<double>, General, RowMajor, Allocator2>& u,
4256  Matrix<complex<double>, General, RowMajor, Allocator3>& v,
4257  LapackInfo& info)
4258  {
4259  int m = A.GetM();
4260  int n = A.GetN();
4261  char jobl('A'), jobr('A');
4262  int lwork = 2 * min(m, n) + max(m, n);
4263  Vector<complex<double> > work(lwork);
4264  Vector<double> rwork(5 * min(m, n));
4265  lambda.Reallocate(min(m, n));
4266  lambda.Zero();
4267  u.Reallocate(m, m);
4268  u.Zero();
4269  v.Reallocate(n, n);
4270  v.Zero();
4271  zgesvd_(&jobl, &jobr, &n, &m, A.GetDataVoid(), &n, lambda.GetData(),
4272  v.GetDataVoid(), &n, u.GetDataVoid(), &m, work.GetDataVoid(),
4273  &lwork, rwork.GetData(), &info.GetInfoRef());
4274 
4275 #ifdef SELDON_LAPACK_CHECK_INFO
4276  if (info.GetInfo() != 0)
4277  throw LapackError(info.GetInfo(), "GetSVD",
4278  "Failed to find singular value decomposition");
4279 #endif
4280 
4281  }
4282 
4283 
4284  /* ColMajor */
4285 
4286 
4287  template<class Prop1, class Allocator1, class Allocator2,
4288  class Allocator3, class Allocator4>
4289  void GetSVD(Matrix<float, Prop1, ColMajor, Allocator1>& A,
4290  Vector<float, VectFull, Allocator4>& lambda,
4291  Matrix<float, General, ColMajor, Allocator2>& u,
4292  Matrix<float, General, ColMajor, Allocator3>& v,
4293  LapackInfo& info)
4294  {
4295  int m = A.GetM();
4296  int n = A.GetN();
4297  char jobl('A'), jobr('A');
4298  int lwork = max(3 * min(m, n) + max(m, n), 5 * min(m, n));
4299  Vector<float> work(lwork);
4300  lambda.Reallocate(min(m, n));
4301  lambda.Zero();
4302  u.Reallocate(m, m);
4303  u.Zero();
4304  v.Reallocate(n, n);
4305  v.Zero();
4306  sgesvd_(&jobl, &jobr, &m, &n, A.GetData(), &m, lambda.GetData(),
4307  u.GetData(), &m, v.GetData(), &n, work.GetData(),
4308  &lwork, &info.GetInfoRef());
4309 
4310 #ifdef SELDON_LAPACK_CHECK_INFO
4311  if (info.GetInfo() != 0)
4312  throw LapackError(info.GetInfo(), "GetSVD",
4313  "Failed to find singular value decomposition");
4314 #endif
4315 
4316  }
4317 
4318 
4319  template<class Prop1, class Allocator1, class Allocator2,
4320  class Allocator3, class Allocator4>
4321  void GetSVD(Matrix<complex<float>, Prop1, ColMajor, Allocator1>& A,
4322  Vector<float, VectFull, Allocator4>& lambda,
4323  Matrix<complex<float>, General, ColMajor, Allocator2>& u,
4324  Matrix<complex<float>, General, ColMajor, Allocator3>& v,
4325  LapackInfo& info)
4326  {
4327  int m = A.GetM();
4328  int n = A.GetN();
4329  char jobl('A'), jobr('A');
4330  int lwork = 2 * min(m, n) + max(m, n);
4331  Vector<complex<float> > work(lwork);
4332  Vector<float> rwork(5 * min(m, n));
4333  lambda.Reallocate(min(m, n));
4334  lambda.Zero();
4335  u.Reallocate(m, m);
4336  u.Zero();
4337  v.Reallocate(n, n);
4338  v.Zero();
4339  cgesvd_(&jobl, &jobr, &m, &n, A.GetDataVoid(), &m, lambda.GetData(),
4340  u.GetDataVoid(), &m, v.GetDataVoid(), &n, work.GetDataVoid(),
4341  &lwork, rwork.GetData(), &info.GetInfoRef());
4342 
4343 #ifdef SELDON_LAPACK_CHECK_INFO
4344  if (info.GetInfo() != 0)
4345  throw LapackError(info.GetInfo(), "GetSVD",
4346  "Failed to find singular value decomposition");
4347 #endif
4348 
4349  }
4350 
4351 
4352  template<class Prop1, class Allocator1, class Allocator2,
4353  class Allocator3, class Allocator4>
4354  void GetSVD(Matrix<double, Prop1, ColMajor, Allocator1>& A,
4355  Vector<double, VectFull, Allocator4>& lambda,
4356  Matrix<double, General, ColMajor, Allocator2>& u,
4357  Matrix<double, General, ColMajor, Allocator3>& v,
4358  LapackInfo& info)
4359  {
4360  int m = A.GetM();
4361  int n = A.GetN();
4362  char jobl('A'), jobr('A');
4363  int lwork = 10 * max(m, n);
4364  Vector<double> work(lwork);
4365  lambda.Reallocate(min(m, n));
4366  lambda.Zero();
4367  u.Reallocate(m, m);
4368  u.Zero();
4369  v.Reallocate(n, n);
4370  v.Zero();
4371  dgesvd_(&jobl, &jobr, &m, &n, A.GetData(), &m, lambda.GetData(),
4372  u.GetData(), &m, v.GetData(), &n, work.GetData(),
4373  &lwork, &info.GetInfoRef());
4374 
4375 #ifdef SELDON_LAPACK_CHECK_INFO
4376  if (info.GetInfo() != 0)
4377  throw LapackError(info.GetInfo(), "GetSVD",
4378  "Failed to find singular value decomposition");
4379 #endif
4380 
4381  }
4382 
4383 
4384  template<class Prop1, class Allocator1, class Allocator2,
4385  class Allocator3, class Allocator4>
4386  void GetSVD(Matrix<complex<double>, Prop1, ColMajor, Allocator1>& A,
4387  Vector<double, VectFull, Allocator4>& sigma,
4388  Matrix<complex<double>, General, ColMajor, Allocator2>& u,
4389  Matrix<complex<double>, General, ColMajor, Allocator3>& v,
4390  LapackInfo& info)
4391  {
4392  int m = A.GetM();
4393  int n = A.GetN();
4394  char jobl('A'), jobr('A');
4395  int lwork = 2 * min(m, n) + max(m, n);
4396  Vector<complex<double> > work(lwork);
4397  Vector<double> rwork(5 * min(m, n));
4398  sigma.Reallocate(min(m, n));
4399  sigma.Zero();
4400  u.Reallocate(m, m);
4401  u.Zero();
4402  v.Reallocate(n, n);
4403  v.Zero();
4404  zgesvd_(&jobl, &jobr, &m, &n, A.GetDataVoid(), &m, sigma.GetData(),
4405  u.GetDataVoid(), &m, v.GetDataVoid(), &n, work.GetDataVoid(),
4406  &lwork, rwork.GetData(), &info.GetInfoRef());
4407 
4408 #ifdef SELDON_LAPACK_CHECK_INFO
4409  if (info.GetInfo() != 0)
4410  throw LapackError(info.GetInfo(), "GetSVD",
4411  "Failed to find singular value decomposition");
4412 #endif
4413 
4414  }
4415 
4416 
4417  // pseudo inverse
4418  template<class T, class Prop, class Storage, class Allocator>
4419  void GetPseudoInverse(Matrix<T, Prop, Storage, Allocator>& A,
4420  const T& epsilon, LapackInfo& info)
4421  {
4422  int m = A.GetM(), n = A.GetN();
4423  Vector<T, VectFull, Allocator> lambda;
4424  Matrix<T, Prop, Storage, Allocator> U;
4425  Matrix<T, Prop, Storage, Allocator> V;
4426 
4427  GetSVD(A, lambda, U, V);
4428 
4429  A.Reallocate(n, m);
4430  A.Fill(0.0);
4431  // computation of A = V* Sigma^{-1} U*
4432  for (int k = 0; k < min(m, n); k++)
4433  if (abs(lambda(k)) > epsilon)
4434  for (int j = 0; j < m; j++)
4435  U(j, k) /= lambda(k);
4436 
4437  U.Resize(m, n);
4438  for (int k = m; k < n; k++)
4439  for (int j = 0; j < m; j++)
4440  U(j, k) = 0.0;
4441 
4442  MltAdd(1.0, SeldonTrans, V, SeldonTrans, U, 0.0, A);
4443  }
4444 
4445 
4446  // pseudo inverse
4447  template<class T, class Prop, class Storage, class Allocator>
4448  void GetPseudoInverse(Matrix<complex<T>, Prop, Storage, Allocator>& A,
4449  const T& epsilon, LapackInfo& info)
4450  {
4451  int m = A.GetM(), n = A.GetN();
4452  Vector<T> lambda;
4453  Matrix<complex<T>, Prop, Storage, Allocator> U;
4454  Matrix<complex<T>, Prop, Storage, Allocator> V;
4455 
4456  GetSVD(A, lambda, U, V, info);
4457 
4458  complex<T> one(1.0, 0.0), zero(0.0, 0.0);
4459 
4460  A.Reallocate(n, m);
4461  A.Fill(zero);
4462  // computation of A = V* Sigma^{-1} U*
4463  for (int k = 0; k < min(m, n); k++)
4464  if (abs(lambda(k)) > epsilon)
4465  for (int j = 0; j < m; j++)
4466  U(j, k) /= lambda(k);
4467 
4468  U.Resize(m, n);
4469  for (int k = m; k < n; k++)
4470  for (int j = 0; j < m; j++)
4471  U(j, k) = zero;
4472 
4473  MltAdd(one, SeldonConjTrans, V, SeldonConjTrans, U, zero, A);
4474  }
4475 
4476 
4477  // SINGULAR VALUE DECOMPOSITION //
4479 
4480 
4482  // RESOLUTION SYLVESTER EQUATION //
4483 
4484 
4486 
4493  template<class Alloc>
4494  void GetHessenberg(Matrix<complex<double>, General, ColMajor, Alloc>& A,
4495  Matrix<complex<double>, General, ColMajor, Alloc>& Q,
4496  LapackInfo& info)
4497  {
4498 #ifdef SELDON_CHECK_DIMENSIONS
4499  if (A.GetM() != A.GetN())
4500  throw WrongDim("GetHessenberg", "Matrix must be squared");
4501 #endif
4502 
4503  int n = A.GetM();
4504  int ilo = 1, ihi = n;
4505  Vector<complex<double>, VectFull, Alloc> tau(n-1);
4506  int lwork = n;
4507  Vector<complex<double>, VectFull, Alloc> work(lwork);
4508  zgehrd_(&n, &ilo, &ihi, A.GetDataVoid(), &n, tau.GetDataVoid(),
4509  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4510 
4511 #ifdef SELDON_LAPACK_CHECK_INFO
4512  if (info.GetInfo() != 0)
4513  throw LapackError(info.GetInfo(), "GetHessenberg",
4514  "Failed to reduce A to Hessenberg form");
4515 #endif
4516 
4517  // generating Q
4518  Q = A;
4519  complex<double> zero(0, 0);
4520  for (int i = 0; i < n; i++)
4521  for (int j = 0; j < i-1; j++)
4522  A(i, j) = zero;
4523 
4524  zunghr_(&n, &ilo, &ihi, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4525  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4526 
4527 #ifdef SELDON_LAPACK_CHECK_INFO
4528  if (info.GetInfo() != 0)
4529  throw LapackError(info.GetInfo(), "GetHessenberg",
4530  "Failed to generate unitary matrix Q");
4531 #endif
4532  }
4533 
4534 
4536 
4545  template<class Alloc>
4546  void GetHessenberg(Matrix<complex<double>, General, ColMajor, Alloc>& A,
4547  Matrix<complex<double>, General, ColMajor, Alloc>& B,
4548  Matrix<complex<double>, General, ColMajor, Alloc>& Q,
4549  Matrix<complex<double>, General, ColMajor, Alloc>& Z,
4550  LapackInfo& info)
4551  {
4552 #ifdef SELDON_CHECK_DIMENSIONS
4553  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
4554  throw WrongDim("GetHessenberg",
4555  "Matrix A and B must be squared");
4556 
4557  if (A.GetM() != B.GetM())
4558  throw WrongDim("GetHessenberg",
4559  "Matrix A and B must have the same size");
4560 #endif
4561 
4562  char compq('V'), compz('I');
4563  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4 * n;
4564  Vector<complex<double> > tau(n);
4565  Vector<complex<double> > work(lwork);
4566  zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4567  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4568 
4569 #ifdef SELDON_LAPACK_CHECK_INFO
4570  if (info.GetInfo() != 0)
4571  throw LapackError(info.GetInfo(), "GetHessenberg",
4572  "Failed to compute QR factorisation of B");
4573 #endif
4574 
4575  Q = B;
4576  zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4577  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4578 
4579 #ifdef SELDON_LAPACK_CHECK_INFO
4580  if (info.GetInfo() != 0)
4581  throw LapackError(info.GetInfo(), "GetHessenberg",
4582  "Failed to generate unitary matrix Q");
4583 #endif
4584 
4585  char side('L'), trans('C');
4586  zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4587  A.GetDataVoid(), &n, work.GetData(), &lwork, &info.GetInfoRef());
4588 
4589  for (int i = 0; i < n; i++)
4590  for (int j = 0; j < i; j++)
4591  B(i, j) = 0;
4592 
4593  Z.Reallocate(n, n);
4594  zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4595  B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(),
4596  &n, &info.GetInfoRef());
4597 
4598 #ifdef SELDON_LAPACK_CHECK_INFO
4599  if (info.GetInfo() != 0)
4600  throw LapackError(info.GetInfo(), "GetHessenberg",
4601  "Failed to generate unitary matrix Z");
4602 #endif
4603 
4604  }
4605 
4606 
4608 
4616  template<class Alloc>
4617  void GetQZ(Matrix<complex<double>, General, ColMajor, Alloc>& A,
4618  Matrix<complex<double>, General, ColMajor, Alloc>& B,
4619  Matrix<complex<double>, General, ColMajor, Alloc>& Q,
4620  Matrix<complex<double>, General, ColMajor, Alloc>& Z,
4621  LapackInfo& info)
4622  {
4623 #ifdef SELDON_CHECK_DIMENSIONS
4624  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
4625  throw WrongDim("GetQZ",
4626  "Matrix A and B must be squared");
4627 
4628  if (A.GetM() != B.GetM())
4629  throw WrongDim("GetQZ",
4630  "Matrix A and B must have the same size");
4631 #endif
4632 
4633  char compq('V'), compz('I');
4634  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4*n;
4635  Vector<complex<double> > tau(n);
4636  Vector<complex<double> > work(lwork);
4637  zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4638  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4639 
4640 #ifdef SELDON_LAPACK_CHECK_INFO
4641  if (info.GetInfo() != 0)
4642  throw LapackError(info.GetInfo(), "GetQZ",
4643  "Failed to compute QR factorisation of B");
4644 #endif
4645 
4646  Q = B;
4647  zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4648  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4649 
4650 #ifdef SELDON_LAPACK_CHECK_INFO
4651  if (info.GetInfo() != 0)
4652  throw LapackError(info.GetInfo(), "GetQZ",
4653  "Failed to generate unitary matrix Q");
4654 #endif
4655 
4656  char side('L'), trans('C');
4657  zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4658  A.GetDataVoid(), &n, work.GetData(), &lwork, &info.GetInfoRef());
4659 
4660  for (int i = 0; i < n; i++)
4661  for (int j = 0; j < i; j++)
4662  B(i,j) = 0;
4663 
4664  zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4665  B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(),
4666  &n, &info.GetInfoRef());
4667 
4668  char job('S');
4669  compq = 'V';
4670  compz = 'V';
4671  Vector<complex<double> > alpha(n), beta(n);
4672  Vector<double> rwork(lwork);
4673  zhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4674  B.GetDataVoid(), &n, alpha.GetDataVoid(), beta.GetDataVoid(),
4675  Q.GetDataVoid(), &n, Z.GetDataVoid(), &n, work.GetDataVoid(),
4676  &lwork, rwork.GetData(), &info.GetInfoRef());
4677 
4678 #ifdef SELDON_LAPACK_CHECK_INFO
4679  if (info.GetInfo() != 0)
4680  throw LapackError(info.GetInfo(), "GetQZ",
4681  "Failed to generate qz factorisation");
4682 #endif
4683  }
4684 
4685 
4686  template<class Alloc>
4687  void GetHessenberg(Matrix<complex<double>, General, RowMajor, Alloc>& A,
4688  Matrix<complex<double>, General, RowMajor, Alloc>& Q,
4689  LapackInfo& info)
4690  {
4691 #ifdef SELDON_CHECK_DIMENSIONS
4692  if (A.GetM() != A.GetN())
4693  throw WrongDim("GetHessenberg", "Matrix must be squared");
4694 #endif
4695 
4696  int n = A.GetM();
4697  int ilo = 1, ihi = n;
4698  Vector<complex<double>, VectFull, Alloc> tau(n-1);
4699  int lwork = n;
4700  Vector<complex<double>, VectFull, Alloc> work(lwork);
4701  Transpose(A);
4702  zgehrd_(&n, &ilo, &ihi, A.GetDataVoid(), &n, tau.GetDataVoid(),
4703  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4704 
4705 #ifdef SELDON_LAPACK_CHECK_INFO
4706  if (info.GetInfo() != 0)
4707  throw LapackError(info.GetInfo(), "GetHessenberg",
4708  "Failed to reduce A to Hessenberg form");
4709 #endif
4710 
4711  // generating Q
4712  Q = A;
4713  zunghr_(&n, &ilo, &ihi, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4714  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4715 
4716 #ifdef SELDON_LAPACK_CHECK_INFO
4717  if (info.GetInfo() != 0)
4718  throw LapackError(info.GetInfo(), "GetHessenberg",
4719  "Failed to generate unitary matrix Q");
4720 #endif
4721 
4722  Transpose(Q);
4723  Transpose(A);
4724  complex<double> zero(0, 0);
4725  for (int i = 0; i < n; i++)
4726  for (int j = 0; j < i-1; j++)
4727  A(i, j) = zero;
4728  }
4729 
4730 
4731  template<class Alloc>
4732  void GetHessenberg(Matrix<complex<double>, General, RowMajor, Alloc>& A,
4733  Matrix<complex<double>, General, RowMajor, Alloc>& B,
4734  Matrix<complex<double>, General, RowMajor, Alloc>& Q,
4735  Matrix<complex<double>, General, RowMajor, Alloc>& Z,
4736  LapackInfo& info)
4737  {
4738 #ifdef SELDON_CHECK_DIMENSIONS
4739  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
4740  throw WrongDim("GetHessenberg",
4741  "Matrix A and B must be squared");
4742 
4743  if (A.GetM() != B.GetM())
4744  throw WrongDim("GetHessenberg",
4745  "Matrix A and B must have the same size");
4746 #endif
4747 
4748  char compq('V'), compz('I');
4749  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4 * n;
4750  Transpose(A); Transpose(B);
4751  Vector<complex<double> > tau(n);
4752  Vector<complex<double> > work(lwork);
4753  zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4754  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4755 
4756 #ifdef SELDON_LAPACK_CHECK_INFO
4757  if (info.GetInfo() != 0)
4758  throw LapackError(info.GetInfo(), "GetHessenberg",
4759  "Failed to compute QR factorisation of B");
4760 #endif
4761 
4762  Q = B;
4763  zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4764  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4765 
4766 #ifdef SELDON_LAPACK_CHECK_INFO
4767  if (info.GetInfo() != 0)
4768  throw LapackError(info.GetInfo(), "GetHessenberg",
4769  "Failed to generate unitary matrix Q");
4770 #endif
4771 
4772  char side('L'), trans('C');
4773  zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4774  A.GetDataVoid(), &n, work.GetData(), &lwork, &info.GetInfoRef());
4775 
4776  for (int i = 0; i < n; i++)
4777  for (int j = 0; j < i; j++)
4778  B(j, i) = 0;
4779 
4780  Z.Reallocate(n, n);
4781  zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4782  B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(),
4783  &n, &info.GetInfoRef());
4784 
4785 #ifdef SELDON_LAPACK_CHECK_INFO
4786  if (info.GetInfo() != 0)
4787  throw LapackError(info.GetInfo(), "GetHessenberg",
4788  "Failed to generate unitary matrix Z");
4789 #endif
4790 
4791  Transpose(A); Transpose(B);
4792  Transpose(Q); Transpose(Z);
4793  }
4794 
4795 
4796  template<class Alloc>
4797  void GetQZ(Matrix<complex<double>, General, RowMajor, Alloc>& A,
4798  Matrix<complex<double>, General, RowMajor, Alloc>& B,
4799  Matrix<complex<double>, General, RowMajor, Alloc>& Q,
4800  Matrix<complex<double>, General, RowMajor, Alloc>& Z,
4801  LapackInfo& info)
4802  {
4803 #ifdef SELDON_CHECK_DIMENSIONS
4804  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
4805  throw WrongDim("GetQZ",
4806  "Matrix A and B must be squared");
4807 
4808  if (A.GetM() != B.GetM())
4809  throw WrongDim("GetQZ",
4810  "Matrix A and B must have the same size");
4811 #endif
4812 
4813  char compq('V'), compz('I');
4814  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4*n;
4815  Transpose(A); Transpose(B);
4816  Vector<complex<double> > tau(n);
4817  Vector<complex<double> > work(lwork);
4818  zgeqrf_(&n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4819  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4820 
4821 #ifdef SELDON_LAPACK_CHECK_INFO
4822  if (info.GetInfo() != 0)
4823  throw LapackError(info.GetInfo(), "GetQZ",
4824  "Failed to compute QR factorisation of B");
4825 #endif
4826 
4827  Q = B;
4828  zungqr_(&n, &n, &n, Q.GetDataVoid(), &n, tau.GetDataVoid(),
4829  work.GetDataVoid(), &lwork, &info.GetInfoRef());
4830 
4831 #ifdef SELDON_LAPACK_CHECK_INFO
4832  if (info.GetInfo() != 0)
4833  throw LapackError(info.GetInfo(), "GetQZ",
4834  "Failed to generate unitary matrix Q");
4835 #endif
4836 
4837  char side('L'), trans('C');
4838  zunmqr_(&side, &trans, &n, &n, &n, B.GetDataVoid(), &n, tau.GetDataVoid(),
4839  A.GetDataVoid(), &n, work.GetData(), &lwork, &info.GetInfoRef());
4840 
4841  for (int i = 0; i < n; i++)
4842  for (int j = 0; j < i; j++)
4843  B(j,i) = 0;
4844 
4845  Z.Reallocate(n, n);
4846  zgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4847  B.GetDataVoid(), &n, Q.GetDataVoid(), &n, Z.GetDataVoid(),
4848  &n, &info.GetInfoRef());
4849 
4850  char job('S');
4851  compq = 'V';
4852  compz = 'V';
4853  Vector<complex<double> > alpha(n), beta(n);
4854  Vector<double> rwork(lwork);
4855  zhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetDataVoid(), &n,
4856  B.GetDataVoid(), &n, alpha.GetDataVoid(), beta.GetDataVoid(),
4857  Q.GetDataVoid(), &n, Z.GetDataVoid(), &n, work.GetDataVoid(),
4858  &lwork, rwork.GetData(), &info.GetInfoRef());
4859 
4860 #ifdef SELDON_LAPACK_CHECK_INFO
4861  if (info.GetInfo() != 0)
4862  throw LapackError(info.GetInfo(), "GetQZ",
4863  "Failed to generate qz factorisation");
4864 #endif
4865 
4866  Transpose(A); Transpose(B);
4867  Transpose(Q); Transpose(Z);
4868  }
4869 
4870 
4872  template<class T, class Prop, class Storage, class Allocator, class Vector1>
4873  void SolveHessenberg(Matrix<T, Prop, Storage, Allocator>& A, Vector1& B)
4874  {
4875  int n = A.GetM();
4876  T tmp, pivot, invDiag;
4877  // loop over rows
4878  for (int i = 0; i < n-1; i++)
4879  {
4880  // pivoting
4881  if (abs(A(i+1, i)) > abs(A(i, i)))
4882  {
4883  // swapping rows
4884  for (int j = i; j < n; j++)
4885  {
4886  tmp = A(i, j);
4887  A(i, j) = A(i+1, j);
4888  A(i+1, j) = tmp;
4889  }
4890 
4891  tmp = B(i);
4892  B(i) = B(i+1);
4893  B(i+1) = tmp;
4894  }
4895 
4896  // performing elimination of A(i+1, i)
4897  invDiag = 1.0/A(i, i);
4898  pivot = A(i+1, i)*invDiag;
4899  A(i, i) = invDiag;
4900  A(i+1, i) = 0;
4901  for (int j = i+1; j < n; j++)
4902  A(i+1, j) -= pivot*A(i, j);
4903 
4904  B(i+1) -= pivot*B(i);
4905  }
4906 
4907  // inverting last element
4908  A(n-1, n-1) = 1.0/A(n-1, n-1);
4909 
4910  // then solving triangular system
4911  for (int i = n-1; i >= 0; i--)
4912  {
4913  tmp = B(i);
4914  for (int j = i+1; j < n; j++)
4915  tmp -= A(i, j)*B(j);
4916 
4917  B(i) = tmp*A(i, i);
4918  }
4919  }
4920 
4921 
4924  template<class T, class Prop, class Storage, class Allocator, class Vector1>
4925  void SolveHessenbergTwo(Matrix<T, Prop, Storage, Allocator>& A, Vector1& B)
4926  {
4927  int n = A.GetM();
4928  T tmp, pivot, invDiag, one, zero;
4929  SetComplexZero(zero);
4930  SetComplexOne(one);
4931  typename ClassComplexType<T>::Treal a1, a2, a3;
4932  // loop over rows
4933  for (int i = 0; i < n-2; i++)
4934  {
4935  // pivoting
4936  a1 = abs(A(i, i));
4937  a2 = abs(A(i+1, i));
4938  a3 = abs(A(i+2, i));
4939  if (a2 > max(a1, a3))
4940  {
4941  // swapping rows i and i+1
4942  for (int j = i; j < n; j++)
4943  {
4944  tmp = A(i, j);
4945  A(i, j) = A(i+1, j);
4946  A(i+1, j) = tmp;
4947  }
4948 
4949  tmp = B(i);
4950  B(i) = B(i+1);
4951  B(i+1) = tmp;
4952  }
4953  else if (a3 > max(a1, a2))
4954  {
4955  // swapping rows i+1 and i+2
4956  for (int j = i; j < n; j++)
4957  {
4958  tmp = A(i, j);
4959  A(i, j) = A(i+2, j);
4960  A(i+2, j) = tmp;
4961  }
4962 
4963  tmp = B(i);
4964  B(i) = B(i+2);
4965  B(i+2) = tmp;
4966  }
4967 
4968  // performing elimination of A(i+1, i)
4969  invDiag = 1.0/A(i, i);
4970  pivot = A(i+1, i)*invDiag;
4971  A(i, i) = invDiag;
4972  A(i+1, i) = zero;
4973  for (int j = i+1; j < n; j++)
4974  A(i+1, j) -= pivot*A(i, j);
4975 
4976  B(i+1) -= pivot*B(i);
4977 
4978  // then elimination of A(i+2, i)
4979  pivot = A(i+2, i)*invDiag;
4980  A(i+2, i) = zero;
4981  for (int j = i+1; j < n; j++)
4982  A(i+2, j) -= pivot*A(i, j);
4983 
4984  B(i+2) -= pivot*B(i);
4985  }
4986 
4987  // elimination of A(n, n-1)
4988  if (abs(A(n-1, n-2)) > abs(A(n-2, n-2)))
4989  {
4990  for (int j = n-2; j < n; j++)
4991  {
4992  tmp = A(n-2, j);
4993  A(n-2, j) = A(n-1, j);
4994  A(n-1, j) = tmp;
4995  }
4996 
4997  tmp = B(n-2);
4998  B(n-2) = B(n-1);
4999  B(n-1) = tmp;
5000  }
5001 
5002  invDiag = one/A(n-2, n-2);
5003  pivot = A(n-1, n-2)*invDiag;
5004  A(n-2, n-2) = invDiag;
5005  A(n-1, n-2) = 0;
5006  A(n-1, n-1) -= pivot*A(n-2, n-1);
5007  B(n-1) -= pivot*B(n-2);
5008 
5009  // inverting last element
5010  A(n-1, n-1) = one/A(n-1, n-1);
5011 
5012  // then solving triangular system
5013  for (int i = n-1; i >= 0; i--)
5014  {
5015  tmp = B(i);
5016  for (int j = i+1; j < n; j++)
5017  tmp -= A(i, j)*B(j);
5018 
5019  B(i) = tmp*A(i, i);
5020  }
5021  }
5022 
5023 
5026  template<class Prop, class Storage, class Allocator>
5027  void SolveSylvester(Matrix<complex<double>, Prop, Storage, Allocator>& A,
5028  Matrix<complex<double>, Prop, Storage, Allocator>& B,
5029  Matrix<complex<double>, Prop, Storage, Allocator>& C,
5030  Matrix<complex<double>, Prop, Storage, Allocator>& D,
5031  Matrix<complex<double>, Prop, Storage, Allocator>& E)
5032  {
5033  complex<double> one(1), zero(0);
5034  int n = A.GetM();
5035 
5036  if (n <= 0)
5037  return;
5038 
5039  if (n == 1)
5040  {
5041  E(0, 0) /= A(0, 0) * conj(B(0, 0)) + C(0, 0) * conj(D(0, 0));
5042  return;
5043  }
5044 
5045  Matrix<complex<double>, Prop, Storage, Allocator> Q1(n, n), Q2(n, n),
5046  Z1(n, n), Z2(n, n);
5047  Matrix<complex<double>, Prop, Storage, Allocator> Y(n, n), F(n, n);
5048 
5049  GetHessenberg(A, C, Q1, Z1);
5050  GetQZ(D, B, Q2, Z2);
5051 
5052  Y.Zero();
5053  MltAdd(one, SeldonConjTrans, Q1, SeldonNoTrans, E, zero, Y);
5054  MltAdd(one, SeldonNoTrans, Y, SeldonNoTrans, Q2, zero, F);
5055 
5056  // Q1 = A y_j, Q2 = C y_j
5057  Vector<complex<double> > Yvec(n);
5058  Q1.Zero(); Q2.Zero(); E.Zero();
5059  complex<double> coef_b, coef_d;
5060  for (int k = n-1; k >= 0; k--)
5061  {
5062  // computation of Yvec = F_k
5063  // - \sum_{j=k+1}^n conj(b_{k, j}) A y_j + conj(d_{k, j}) C y_j
5064  for (int j = 0; j < n; j++)
5065  Yvec(j) = F(j, k);
5066 
5067  for (int j = k+1; j < n; j++)
5068  {
5069  coef_b = conj(B(k, j));
5070  coef_d = conj(D(k, j));
5071  for (int i = 0; i < n; i++)
5072  Yvec(i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j);
5073  }
5074 
5075  coef_b = conj(B(k, k)); coef_d = conj(D(k, k));
5076  for (int i = 0; i < n; i++)
5077  for (int j = max(0, i-1); j < n; j++)
5078  E(i, j) = coef_b * A(i, j) + coef_d * C(i, j);
5079 
5080  SolveHessenberg(E, Yvec);
5081  for (int i = 0; i < n; i++)
5082  Y(i, k) = Yvec(i);
5083 
5084  // computation of A y_k and C y_k
5085  for (int i = 0; i < n; i++)
5086  for (int m = max(0, i-1); m < n; m++)
5087  Q1(i, k) += A(i, m)*Yvec(m);
5088 
5089  for (int i = 0; i < n; i++)
5090  for (int m = i; m < n; m++)
5091  Q2(i, k) += C(i, m)*Yvec(m);
5092  }
5093 
5094  MltAdd(one, SeldonNoTrans, Y, SeldonConjTrans, Z2, zero, F);
5095  MltAdd(one, Z1, F, zero, E);
5096  }
5097 
5098 
5099  template<class Alloc>
5100  void GetHessenberg(Matrix<double, General, ColMajor, Alloc>& A,
5101  Matrix<double, General, ColMajor, Alloc>& Q,
5102  LapackInfo& info)
5103  {
5104 #ifdef SELDON_CHECK_DIMENSIONS
5105  if (A.GetM() != A.GetN())
5106  throw WrongDim("GetHessenberg", "Matrix must be squared");
5107 #endif
5108 
5109  int n = A.GetM();
5110  int ilo = 1, ihi = n;
5111  Vector<double, VectFull, Alloc> tau(n-1);
5112  int lwork = n;
5113  Vector<double, VectFull, Alloc> work(lwork);
5114  dgehrd_(&n, &ilo, &ihi, A.GetData(), &n, tau.GetData(),
5115  work.GetData(), &lwork, &info.GetInfoRef());
5116 
5117 #ifdef SELDON_LAPACK_CHECK_INFO
5118  if (info.GetInfo() != 0)
5119  throw LapackError(info.GetInfo(), "GetHessenberg",
5120  "Failed to reduce A to Hessenberg form");
5121 #endif
5122 
5123  // generating Q
5124  Q = A;
5125  double zero(0);
5126  for (int i = 0; i < n; i++)
5127  for (int j = 0; j < i-1; j++)
5128  A(i, j) = zero;
5129 
5130  dorghr_(&n, &ilo, &ihi, Q.GetData(), &n, tau.GetData(),
5131  work.GetData(), &lwork, &info.GetInfoRef());
5132 
5133 #ifdef SELDON_LAPACK_CHECK_INFO
5134  if (info.GetInfo() != 0)
5135  throw LapackError(info.GetInfo(), "GetHessenberg",
5136  "Failed to generate unitary matrix Q");
5137 #endif
5138  }
5139 
5140 
5142 
5151  template<class Alloc>
5152  void GetHessenberg(Matrix<double, General, ColMajor, Alloc>& A,
5153  Matrix<double, General, ColMajor, Alloc>& B,
5154  Matrix<double, General, ColMajor, Alloc>& Q,
5155  Matrix<double, General, ColMajor, Alloc>& Z,
5156  LapackInfo& info)
5157  {
5158 #ifdef SELDON_CHECK_DIMENSIONS
5159  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
5160  throw WrongDim("GetHessenberg",
5161  "Matrix A and B must be squared");
5162 
5163  if (A.GetM() != B.GetM())
5164  throw WrongDim("GetHessenberg",
5165  "Matrix A and B must have the same size");
5166 #endif
5167 
5168  char compq('V'), compz('I');
5169  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4 * n;
5170  Vector<double> tau(n);
5171  Vector<double> work(lwork);
5172  dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(),
5173  work.GetData(), &lwork, &info.GetInfoRef());
5174 
5175 #ifdef SELDON_LAPACK_CHECK_INFO
5176  if (info.GetInfo() != 0)
5177  throw LapackError(info.GetInfo(), "GetHessenberg",
5178  "Failed to compute QR factorisation of B");
5179 #endif
5180 
5181  Q = B;
5182  dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(),
5183  work.GetData(), &lwork, &info.GetInfoRef());
5184 
5185  char side('L'), trans('T');
5186  dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(),
5187  A.GetData(), &n, work.GetData(), &lwork, &info.GetInfoRef());
5188 
5189 #ifdef SELDON_LAPACK_CHECK_INFO
5190  if (info.GetInfo() != 0)
5191  throw LapackError(info.GetInfo(), "GetHessenberg",
5192  "Failed to generate unitary matrix Q");
5193 #endif
5194 
5195  for (int i = 0; i < n; i++)
5196  for (int j = 0; j < i; j++)
5197  B(i, j) = 0;
5198 
5199  Z.Reallocate(n, n);
5200  dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5201  B.GetData(), &n, Q.GetData(), &n, Z.GetData(),
5202  &n, &info.GetInfoRef());
5203 
5204 #ifdef SELDON_LAPACK_CHECK_INFO
5205  if (info.GetInfo() != 0)
5206  throw LapackError(info.GetInfo(), "GetHessenberg",
5207  "Failed to generate unitary matrix Z");
5208 #endif
5209 
5210  }
5211 
5212 
5214 
5222  template<class Alloc>
5223  void GetQZ(Matrix<double, General, ColMajor, Alloc>& A,
5224  Matrix<double, General, ColMajor, Alloc>& B,
5225  Matrix<double, General, ColMajor, Alloc>& Q,
5226  Matrix<double, General, ColMajor, Alloc>& Z,
5227  LapackInfo& info)
5228  {
5229 #ifdef SELDON_CHECK_DIMENSIONS
5230  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
5231  throw WrongDim("GetQZ",
5232  "Matrix A and B must be squared");
5233 
5234  if (A.GetM() != B.GetM())
5235  throw WrongDim("GetQZ",
5236  "Matrix A and B must have the same size");
5237 #endif
5238 
5239  char compq('V'), compz('I');
5240  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4*n;
5241  Vector<double> tau(n);
5242  Vector<double> work(lwork);
5243  dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(),
5244  work.GetData(), &lwork, &info.GetInfoRef());
5245 
5246 #ifdef SELDON_LAPACK_CHECK_INFO
5247  if (info.GetInfo() != 0)
5248  throw LapackError(info.GetInfo(), "GetQZ",
5249  "Failed to compute QR factorisation of B");
5250 #endif
5251 
5252  Q = B;
5253  dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(),
5254  work.GetData(), &lwork, &info.GetInfoRef());
5255 
5256 #ifdef SELDON_LAPACK_CHECK_INFO
5257  if (info.GetInfo() != 0)
5258  throw LapackError(info.GetInfo(), "GetQZ",
5259  "Failed to generate unitary matrix Q");
5260 #endif
5261 
5262  char side('L'), trans('T');
5263  dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(),
5264  A.GetData(), &n, work.GetData(), &lwork, &info.GetInfoRef());
5265 
5266  for (int i = 0; i < n; i++)
5267  for (int j = 0; j < i; j++)
5268  B(i,j) = 0;
5269 
5270  Z.Reallocate(n, n);
5271  dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5272  B.GetData(), &n, Q.GetData(), &n, Z.GetData(),
5273  &n, &info.GetInfoRef());
5274 
5275  char job('S');
5276  compq = 'V';
5277  compz = 'V';
5278  Vector<double> alphar(n), alphai(n), beta(n);
5279  Vector<double> rwork(lwork);
5280  dhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5281  B.GetData(), &n, alphar.GetData(), alphai.GetData(), beta.GetData(),
5282  Q.GetData(), &n, Z.GetData(), &n, work.GetData(),
5283  &lwork, &info.GetInfoRef());
5284 
5285 #ifdef SELDON_LAPACK_CHECK_INFO
5286  if (info.GetInfo() != 0)
5287  throw LapackError(info.GetInfo(), "GetQZ",
5288  "Failed to generate qz factorisation");
5289 #endif
5290 
5291  }
5292 
5293 
5294  template<class Alloc>
5295  void GetHessenberg(Matrix<double, General, RowMajor, Alloc>& A,
5296  Matrix<double, General, RowMajor, Alloc>& Q,
5297  LapackInfo& info)
5298  {
5299 #ifdef SELDON_CHECK_DIMENSIONS
5300  if (A.GetM() != A.GetN())
5301  throw WrongDim("GetHessenberg", "Matrix must be squared");
5302 #endif
5303 
5304  int n = A.GetM();
5305  int ilo = 1, ihi = n;
5306  Vector<double, VectFull, Alloc> tau(n-1);
5307  int lwork = n;
5308  Vector<double, VectFull, Alloc> work(lwork);
5309  Transpose(A);
5310  dgehrd_(&n, &ilo, &ihi, A.GetData(), &n, tau.GetData(),
5311  work.GetData(), &lwork, &info.GetInfoRef());
5312 
5313 #ifdef SELDON_LAPACK_CHECK_INFO
5314  if (info.GetInfo() != 0)
5315  throw LapackError(info.GetInfo(), "GetHessenberg",
5316  "Failed to reduce A to Hessenberg form");
5317 #endif
5318 
5319  // generating Q
5320  Q = A;
5321  dorghr_(&n, &ilo, &ihi, Q.GetData(), &n, tau.GetData(),
5322  work.GetData(), &lwork, &info.GetInfoRef());
5323 
5324 #ifdef SELDON_LAPACK_CHECK_INFO
5325  if (info.GetInfo() != 0)
5326  throw LapackError(info.GetInfo(), "GetHessenberg",
5327  "Failed to generate unitary matrix Q");
5328 #endif
5329 
5330  Transpose(A);
5331  Transpose(Q);
5332  double zero(0);
5333  for (int i = 0; i < n; i++)
5334  for (int j = 0; j < i-1; j++)
5335  A(i, j) = zero;
5336  }
5337 
5338 
5339  template<class Alloc>
5340  void GetHessenberg(Matrix<double, General, RowMajor, Alloc>& A,
5341  Matrix<double, General, RowMajor, Alloc>& B,
5342  Matrix<double, General, RowMajor, Alloc>& Q,
5343  Matrix<double, General, RowMajor, Alloc>& Z,
5344  LapackInfo& info)
5345  {
5346 #ifdef SELDON_CHECK_DIMENSIONS
5347  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
5348  throw WrongDim("GetHessenberg",
5349  "Matrix A and B must be squared");
5350 
5351  if (A.GetM() != B.GetM())
5352  throw WrongDim("GetHessenberg",
5353  "Matrix A and B must have the same size");
5354 #endif
5355 
5356  char compq('V'), compz('I');
5357  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4 * n;
5358  Transpose(A); Transpose(B);
5359  Vector<double> tau(n);
5360  Vector<double> work(lwork);
5361  dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(),
5362  work.GetData(), &lwork, &info.GetInfoRef());
5363 
5364 #ifdef SELDON_LAPACK_CHECK_INFO
5365  if (info.GetInfo() != 0)
5366  throw LapackError(info.GetInfo(), "GetHessenberg",
5367  "Failed to compute QR factorisation of B");
5368 #endif
5369 
5370  Q = B;
5371  dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(),
5372  work.GetData(), &lwork, &info.GetInfoRef());
5373 
5374  char side('L'), trans('T');
5375  dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(),
5376  A.GetData(), &n, work.GetData(), &lwork, &info.GetInfoRef());
5377 
5378  for (int i = 0; i < n; i++)
5379  for (int j = 0; j < i; j++)
5380  B(j, i) = 0;
5381 
5382  Z.Reallocate(n, n);
5383  dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5384  B.GetData(), &n, Q.GetData(), &n, Z.GetData(),
5385  &n, &info.GetInfoRef());
5386 
5387 #ifdef SELDON_LAPACK_CHECK_INFO
5388  if (info.GetInfo() != 0)
5389  throw LapackError(info.GetInfo(), "GetHessenberg",
5390  "Failed to generate unitary matrix Z");
5391 #endif
5392 
5393  Transpose(A); Transpose(B);
5394  Transpose(Q); Transpose(Z);
5395 
5396  }
5397 
5398 
5399  template<class Alloc>
5400  void GetQZ(Matrix<double, General, RowMajor, Alloc>& A,
5401  Matrix<double, General, RowMajor, Alloc>& B,
5402  Matrix<double, General, RowMajor, Alloc>& Q,
5403  Matrix<double, General, RowMajor>& Z,
5404  LapackInfo& info)
5405  {
5406 #ifdef SELDON_CHECK_DIMENSIONS
5407  if ((A.GetM() != A.GetN()) || (B.GetM() != B.GetN()))
5408  throw WrongDim("GetQZ",
5409  "Matrix A and B must be squared");
5410 
5411  if (A.GetM() != B.GetM())
5412  throw WrongDim("GetQZ",
5413  "Matrix A and B must have the same size");
5414 #endif
5415 
5416  char compq('V'), compz('I');
5417  int n = A.GetM(), ilo = 1, ihi = n, lwork = 4*n;
5418  Transpose(A); Transpose(B);
5419  Vector<double> tau(n);
5420  Vector<double> work(lwork);
5421  dgeqrf_(&n, &n, B.GetData(), &n, tau.GetData(),
5422  work.GetData(), &lwork, &info.GetInfoRef());
5423 
5424 #ifdef SELDON_LAPACK_CHECK_INFO
5425  if (info.GetInfo() != 0)
5426  throw LapackError(info.GetInfo(), "GetQZ",
5427  "Failed to compute QR factorisation of B");
5428 #endif
5429 
5430  Q = B;
5431  dorgqr_(&n, &n, &n, Q.GetData(), &n, tau.GetData(),
5432  work.GetData(), &lwork, &info.GetInfoRef());
5433 
5434  char side('L'), trans('T');
5435  dormqr_(&side, &trans, &n, &n, &n, B.GetData(), &n, tau.GetData(),
5436  A.GetData(), &n, work.GetData(), &lwork, &info.GetInfoRef());
5437 
5438  for (int i = 0; i < n; i++)
5439  for (int j = 0; j < i; j++)
5440  B(j, i) = 0;
5441 
5442  Z.Reallocate(n, n);
5443  dgghrd_(&compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5444  B.GetData(), &n, Q.GetData(), &n, Z.GetData(),
5445  &n, &info.GetInfoRef());
5446 
5447  char job('S');
5448  compq = 'V';
5449  compz = 'V';
5450  Vector<double> alphar(n), alphai(n), beta(n);
5451  Vector<double> rwork(lwork);
5452  dhgeqz_(&job, &compq, &compz, &n, &ilo, &ihi, A.GetData(), &n,
5453  B.GetData(), &n, alphar.GetData(), alphai.GetData(), beta.GetData(),
5454  Q.GetData(), &n, Z.GetData(), &n, work.GetData(),
5455  &lwork, &info.GetInfoRef());
5456 
5457 #ifdef SELDON_LAPACK_CHECK_INFO
5458  if (info.GetInfo() != 0)
5459  throw LapackError(info.GetInfo(), "GetQZ",
5460  "Failed to generate qz factorisation");
5461 #endif
5462 
5463  Transpose(A); Transpose(B);
5464  Transpose(Q); Transpose(Z);
5465  }
5466 
5467 
5470  template<class Prop, class Storage, class Allocator>
5471  void SolveSylvester(Matrix<double, Prop, Storage, Allocator>& A,
5472  Matrix<double, Prop, Storage, Allocator>& B,
5473  Matrix<double, Prop, Storage, Allocator>& C,
5474  Matrix<double, Prop, Storage, Allocator>& D,
5475  Matrix<double, Prop, Storage, Allocator>& E)
5476  {
5477  double one(1), zero(0);
5478  int n = A.GetM();
5479  if (n <= 0)
5480  return;
5481 
5482  if (n == 1)
5483  {
5484  E(0,0) /= A(0, 0) * B(0, 0) + C(0, 0) * D(0, 0);
5485  return;
5486  }
5487 
5488  Matrix<double, Prop, Storage, Allocator> Q1(n, n), Z1(n, n);
5489  Matrix<double, Prop, Storage, Allocator> Q2(n, n), Z2(n, n);
5490  Matrix<double, Prop, Storage, Allocator> Y(n, n), F(n, n);
5491 
5492  GetHessenberg(A, C, Q1, Z1);
5493  GetQZ(D, B, Q2, Z2);
5494 
5495  Y.Zero();
5496  MltAdd(one, SeldonTrans, Q1, SeldonNoTrans, E, zero, Y);
5497  MltAdd(one, SeldonNoTrans, Y, SeldonNoTrans, Q2, zero, F);
5498  Y.Zero();
5499 
5500  // Q1 = A y_j, Q2 = C y_j
5501  Vector<double> Yvec(n);
5502  Q1.Zero(); Q2.Zero(); E.Zero();
5503  double coef_b, coef_d, coef_bm1, coef_dm1, coef_b2, coef_d2, coef_d3;
5504  Vector<double> Yr(2*n);
5505  Matrix<double, Prop, Storage, Allocator> Er(2*n, 2*n);
5506  Yr.Zero(); Er.Zero();
5507  int k = n-1;
5508  while (k >= 0)
5509  {
5510  if ((k==0) || (D(k, k-1) == 0))
5511  {
5512  // 1x1 block : same case as for complex matrix
5513 
5514  // computation of Yvec = F_k
5515  // - \sum_{j=k+1}^n b_{k, j} A y_j + d_{k, j} C y_j
5516  for (int j = 0; j < n; j++)
5517  Yvec(j) = F(j, k);
5518 
5519  for (int j = k+1; j < n; j++)
5520  {
5521  coef_b = B(k, j);
5522  coef_d = D(k, j);
5523  for (int i = 0; i < n; i++)
5524  Yvec(i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j);
5525  }
5526 
5527  coef_b = B(k, k); coef_d = D(k, k);
5528  for (int i = 0; i < n; i++)
5529  for (int j = max(0, i-1); j < n; j++)
5530  E(i, j) = coef_b * A(i, j) + coef_d * C(i, j);
5531 
5532  SolveHessenberg(E, Yvec);
5533  for (int i = 0; i < n; i++)
5534  Y(i, k) = Yvec(i);
5535 
5536  // computation of A y_k and C y_k
5537  for (int i = 0; i < n; i++)
5538  for (int m = max(0, i-1); m < n; m++)
5539  Q1(i, k) += A(i, m)*Yvec(m);
5540 
5541  for (int i = 0; i < n; i++)
5542  for (int m = i; m < n; m++)
5543  Q2(i, k) += C(i, m)*Yvec(m);
5544 
5545  k--;
5546  }
5547  else
5548  {
5549  // 2x2 block, we have to solve 2n x 2n linear system
5550 
5551  // computation of Yr = (f_k-1, f_k)
5552  for (int j = 0; j < n; j++)
5553  {
5554  Yr(2*j) = F(j, k-1);
5555  Yr(2*j+1) = F(j, k);
5556  }
5557 
5558  for (int j = k+1; j < n; j++)
5559  {
5560  coef_b = B(k-1, j);
5561  coef_d = D(k-1, j);
5562  for (int i = 0; i < n; i++)
5563  Yr(2*i) -= coef_b * Q1(i, j) + coef_d * Q2(i, j);
5564 
5565  coef_b = B(k, j);
5566  coef_d = D(k, j);
5567  for (int i = 0; i < n; i++)
5568  Yr(2*i+1) -= coef_b * Q1(i, j) + coef_d * Q2(i, j);
5569  }
5570 
5571  // computation of linear system
5572  // b_{k-1, k-1} A + d_{k-1, k-1} C b_{k-1,k} A + d_{k-1,k} C
5573  // d_{k, k-1} C b_{k,k} A + d_{k,k} C
5574  coef_bm1 = B(k-1, k-1); coef_dm1 = D(k-1, k-1);
5575  coef_b2 = B(k-1, k); coef_d2 = D(k-1, k); coef_d3 = D(k, k-1);
5576  coef_b = B(k, k); coef_d = D(k, k);
5577  for (int i = 0; i < n; i++)
5578  for (int j = max(0, i-1); j < n; j++)
5579  {
5580  Er(2*i, 2*j) = coef_bm1 * A(i, j) + coef_dm1 * C(i, j);
5581  Er(2*i, 2*j+1) = coef_b2 * A(i, j) + coef_d2 * C(i, j);
5582  Er(2*i+1, 2*j) = coef_d3 * C(i, j);
5583  Er(2*i+1, 2*j+1) = coef_b * A(i, j) + coef_d * C(i, j);
5584  }
5585 
5586  SolveHessenbergTwo(Er, Yr);
5587 
5588  for (int i = 0; i < n; i++)
5589  {
5590  Y(i, k-1) = Yr(2*i);
5591  Y(i, k) = Yr(2*i+1);
5592  }
5593 
5594  // computation of A y_k and C y_k
5595  for (int i = 0; i < n; i++)
5596  for (int m = max(0, i-1); m < n; m++)
5597  Q1(i, k-1) += A(i, m)*Yr(2*m);
5598 
5599  for (int i = 0; i < n; i++)
5600  for (int m = i; m < n; m++)
5601  Q2(i, k-1) += C(i, m)*Yr(2*m);
5602 
5603  for (int i = 0; i < n; i++)
5604  for (int m = max(0, i-1); m < n; m++)
5605  Q1(i, k) += A(i, m)*Yr(2*m+1);
5606 
5607  for (int i = 0; i < n; i++)
5608  for (int m = i; m < n; m++)
5609  Q2(i, k) += C(i, m)*Yr(2*m+1);
5610 
5611  k -= 2;
5612  }
5613  }
5614 
5615  MltAdd(one, SeldonNoTrans, Y, SeldonTrans, Z2, zero, F);
5616  MltAdd(one, Z1, F, zero, E);
5617  }
5618 
5619 
5620  // RESOLUTION SYLVESTER EQUATION //
5622 
5623 
5624 } // end namespace
5625 
5626 #define SELDON_FILE_LAPACK_EIGENVALUES_CXX
5627 #endif
Seldon::TransposeConj
void TransposeConj(const Matrix< T, Prop, Storage, Allocator > &A, Matrix< T, Prop, Storage, Allocator > &B)
Matrix transposition and conjugation.
Definition: Functions_Matrix.cxx:2925
Seldon::LapackInfo
Definition: Errors.hxx:243
Seldon::Vector
Definition: SeldonHeader.hxx:207
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::GetQZ
void GetQZ(Matrix< complex< double >, General, ColMajor, Alloc > &A, Matrix< complex< double >, General, ColMajor, Alloc > &B, Matrix< complex< double >, General, ColMajor, Alloc > &Q, Matrix< complex< double >, General, ColMajor, Alloc > &Z, LapackInfo &info)
Reduces A and B to quasi-triangular matrices.
Definition: Lapack_Eigenvalues.cxx:4617
Seldon::VectFull
Definition: StorageInline.cxx:74
Seldon::Conjugate
void Conjugate(Matrix< T, Prop, Storage, Allocator > &A)
A is replaced by its conjugate.
Definition: Functions_Matrix.cxx:2915
Seldon::General
Definition: Properties.hxx:26
Seldon::SolveHessenberg
void SolveHessenberg(Matrix< T, Prop, Storage, Allocator > &A, Vector1 &B)
Gaussian elimination to solve A X = B with A an Hessenberg matrix.
Definition: Lapack_Eigenvalues.cxx:4873
Seldon::GetHessenberg
void GetHessenberg(Matrix< complex< double >, General, ColMajor, Alloc > &A, Matrix< complex< double >, General, ColMajor, Alloc > &Q, LapackInfo &info)
Reduces A to their Hessenberg form.
Definition: Lapack_Eigenvalues.cxx:4494
Seldon::Transpose
void Transpose(Matrix< T, Prop, Storage, Allocator > &A)
Matrix transposition.
Definition: Functions_Matrix.cxx:2699
Seldon::SolveSylvester
void SolveSylvester(Matrix< complex< double >, Prop, Storage, Allocator > &A, Matrix< complex< double >, Prop, Storage, Allocator > &B, Matrix< complex< double >, Prop, Storage, Allocator > &C, Matrix< complex< double >, Prop, Storage, Allocator > &D, Matrix< complex< double >, Prop, Storage, Allocator > &E)
Solves A X B^H + C X D^H = E, E is overwritten with result X. A, B, C and D are modified.
Definition: Lapack_Eigenvalues.cxx:5027
Seldon::SolveHessenbergTwo
void SolveHessenbergTwo(Matrix< T, Prop, Storage, Allocator > &A, Vector1 &B)
Gaussian elimination to solve A X = B with A matrix so that a_ij = 0 for i > j+2.
Definition: Lapack_Eigenvalues.cxx:4925
Seldon::WrongDim
Definition: Errors.hxx:102
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::LapackError
Definition: Errors.hxx:225
Seldon::ColMajor
Definition: Storage.hxx:33