CVM Class Library  8.1
This C++ class library encapsulates concepts of vector and different matrices including square, band, symmetric and hermitian ones in Euclidean space of real and complex numbers.
 All Classes Files Functions Variables Typedefs Friends Macros Pages
cvector.cpp
Go to the documentation of this file.
1 // CVM Class Library
2 // http://cvmlib.com
3 //
4 // Copyright Sergei Nikolaev 1992-2014
5 // Distributed under the Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8 
9 #include "cvm.h"
10 #include "blas.h"
11 
13 
14 template<>
15 CVM_API std::complex<float> __dotu<std::complex<float> > (const std::complex<float>* mpd, tint mnSize, tint mnIncr, const std::complex<float>* pd, tint nIncr)
16 {
17  std::complex<float> cRes;
18  VCDOTU (&cRes, &mnSize, mpd, &mnIncr, pd, &nIncr);
19  return cRes;
20 }
21 
22 template<>
23 CVM_API std::complex<double> __dotu<std::complex<double> > (const std::complex<double>* mpd, tint mnSize, tint mnIncr, const std::complex<double>* pd, tint nIncr)
24 {
25  std::complex<double> cRes;
26  VZDOTU (&cRes, &mnSize, mpd, &mnIncr, pd, &nIncr);
27  return cRes;
28 }
29 
30 template<>
31 CVM_API std::complex<float> __dotc<std::complex<float> > (const std::complex<float>* mpd, tint mnSize, tint mnIncr, const std::complex<float>* pd, tint nIncr)
32 {
33  std::complex<float> cRes;
34  VCDOTC (&cRes, &mnSize, mpd, &mnIncr, pd, &nIncr);
35  return cRes;
36 }
37 
38 template<>
39 CVM_API std::complex<double> __dotc<std::complex<double> > (const std::complex<double>* mpd, tint mnSize, tint mnIncr, const std::complex<double>* pd, tint nIncr)
40 {
41  std::complex<double> cRes;
42  VZDOTC (&cRes, &mnSize, mpd, &mnIncr, pd, &nIncr);
43  return cRes;
44 }
45 
47 template<>
48 CVM_API void
50  (bool bLeft,
52  std::complex<float> dAlpha,
54  std::complex<float> dBeta,
56 {
57  CGEMV (bLeft ? Chars::pT() : Chars::pN(),
58 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
59  1,
60 #endif
61  m._pm(), m._pn(), &dAlpha, m._pd(), m._pldm(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
62 }
63 
64 template<>
65 CVM_API void
67  (bool bLeft,
68  const basic_cmatrix<double, std::complex<double> >& m,
69  std::complex<double> dAlpha,
71  std::complex<double> dBeta,
73 {
74  ZGEMV (bLeft ? Chars::pT() : Chars::pN(),
75 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
76  1,
77 #endif
78  m._pm(), m._pn(), &dAlpha, m._pd(), m._pldm(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
79 }
80 
81 template<>
82 CVM_API void
84  (bool bLeft,
85  const basic_scbmatrix<float, std::complex<float> >& m,
86  std::complex<float> dAlpha,
88  std::complex<float> dBeta,
90 {
91  CGBMV (bLeft ? Chars::pT() : Chars::pN(),
92 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
93  1,
94 #endif
95  m._pm(), m._pn(), m._pl(), m._pu(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
96 }
97 
98 template<>
99 CVM_API void
101  (bool bLeft,
102  const basic_scbmatrix<double, std::complex<double> >& m,
103  std::complex<double> dAlpha,
105  std::complex<double> dBeta,
107 {
108  ZGBMV (bLeft ? Chars::pT() : Chars::pN(),
109 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
110  1,
111 #endif
112  m._pm(), m._pn(), m._pl(), m._pu(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
113 }
114 
115 template<>
116 CVM_API void
118  (const basic_schmatrix<float, std::complex<float> >& m,
119  std::complex<float> cAlpha,
121  std::complex<float> cBeta,
123 {
124  // calling with 'L' for left-sided multiplication does not work here
125  CHEMV (Chars::pU(),
126 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
127  1,
128 #endif
129  m._pm(), &cAlpha, m, m._pld(), v, v._pincr(), &cBeta, vRes, vRes._pincr());
130 }
131 
132 template<>
133 CVM_API void
135  (const basic_schmatrix<double, std::complex<double> >& m,
136  std::complex<double> cAlpha,
138  std::complex<double> cBeta,
140 {
141  // calling with 'L' for left-sided multiplication does not work here
142  ZHEMV (Chars::pU(),
143 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
144  1,
145 #endif
146  m._pm(), &cAlpha, m, m._pld(), v, v._pincr(), &cBeta, vRes, vRes._pincr());
147 }
148 
149 template<>
150 CVM_API void
151 __eig<basic_cvector<float, std::complex<float> >, basic_srmatrix<float>, basic_scmatrix<float, std::complex<float> > >
153  const basic_srmatrix<float>& mArg,
155  bool bRightVect) throw (cvmexception)
156 {
157  const bool bEigVect = (mEigVect != nullptr);
158  const tint nM = mArg.msize();
159  tint lWork = -1;
160  tint ilo = 0;
161  tint ihi = 0;
162  tint nOutInfo = 0;
163 
164  _check_ne(CVM_SIZESMISMATCH, vRes.size(), nM);
165  if (nM == 1)
166  {
167  const std::complex<float> one(1.F, 0.F);
168  vRes[CVM0] = std::complex<float>(mArg(CVM0,CVM0), 0.F);
169  if (bEigVect)
170  {
171  mEigVect -> resize (1);
172  (*mEigVect)[CVM0].set(one);
173  }
174  }
175  else
176  {
177  basic_srmatrix<float> mA (mArg);
178  basic_rvector<float> vScale (nM);
179  basic_rvector<float> vTau (_cvm_max<tint>(1, nM - 1));
180  basic_rvector<float> vR (nM);
181  basic_rvector<float> vI (nM);
182 
183  SGEBAL (Chars::pB(),
184 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
185  1,
186 #endif
187  &nM, mA, &nM, &ilo, &ihi, vScale, &nOutInfo);
188  _check_negative(CVM_WRONGMKLARG, nOutInfo);
189 
190  float dWork;
191  SGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, &dWork, &lWork, &nOutInfo);
192  _check_negative(CVM_WRONGMKLARG, nOutInfo);
193  lWork = static_cast<tint>(dWork);
194  basic_rvector<float> vWork (lWork);
195 
196  SGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, vWork, &lWork, &nOutInfo);
197  _check_negative(CVM_WRONGMKLARG, nOutInfo);
198 
199  if (bEigVect)
200  {
201  const tint one(1);
202  tint m = 0;
203  tint lSelect = 0;
204  const tint ldvl = bRightVect ? 1 : nM;
205  const tint ldvr = bRightVect ? nM : 1;
206  basic_srmatrix<float> vl (ldvl);
207  basic_srmatrix<float> vr (ldvr);
208  basic_rvector <float> work (3 * nM);
209  basic_srmatrix<float> mH (mA);
210  const char* pRL = bRightVect ? Chars::pR() : Chars::pL();
211 
212  if (bRightVect)
213  {
214  vr = mA;
215  }
216  else
217  {
218  vl = mA;
219  }
220 
221  lWork = -1;
222  SORGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, &dWork, &lWork, &nOutInfo);
223  _check_negative(CVM_WRONGMKLARG, nOutInfo);
224  lWork = static_cast<tint>(dWork);
225  if (lWork > vWork.size()) vWork.resize(lWork);
226 
227  SORGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, vWork, &lWork, &nOutInfo);
228  _check_negative(CVM_WRONGMKLARG, nOutInfo);
229 
230  lWork = -1;
231 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
232  SHSEQR (Chars::pS(), 1, Chars::pV(), 1,
233 #else
234  SHSEQR (Chars::pS(), Chars::pV(),
235 #endif
236  &nM, &ilo, &ihi, mH, &nM, vR, vI, bRightVect ? vr : vl, &nM, &dWork, &lWork, &nOutInfo);
237  _check_negative(CVM_WRONGMKLARG, nOutInfo);
238  lWork = static_cast<tint>(dWork) * 11;
239  if (lWork > vWork.size()) vWork.resize(lWork);
240 
241 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
242  SHSEQR (Chars::pS(), 1, Chars::pV(), 1,
243 #else
244  SHSEQR (Chars::pS(), Chars::pV(),
245 #endif
246  &nM, &ilo, &ihi, mH, &nM, vR, vI, bRightVect ? vr : vl, &nM, vWork, &lWork, &nOutInfo);
247  _check_negative(CVM_WRONGMKLARG, nOutInfo);
248  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SHSEQR", __FILE__, __LINE__);
249 
250 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
251  STREVC (pRL, 1, Chars::pB(), 1,
252 #else
253  STREVC (pRL, Chars::pB(),
254 #endif
255  &lSelect, &nM, mH, &nM, vl, &ldvl, vr, &ldvr, &nM, &m, work, &nOutInfo);
256  _check_negative(CVM_WRONGMKLARG, nOutInfo);
257 
258  basic_srmatrix<float>& v = bRightVect ? vr : vl;
259  const tint ldv = bRightVect ? ldvr : ldvl;
260 
261 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
262  SGEBAK (Chars::pB(), 1, pRL, 1,
263 #else
264  SGEBAK (Chars::pB(), pRL,
265 #endif
266  &nM, &ilo, &ihi, vScale, &nM, v, &ldv, &nOutInfo);
267  _check_negative(CVM_WRONGMKLARG, nOutInfo);
268 
269  bool bPair = false;
270  m = CVM0;
271  mEigVect -> resize (nM);
272  for (tint i = CVM0; i < nM + CVM0; i++)
273  {
274  if (_abs(vI(i)) > basic_cvmMachMin<float>())
275  {
276  (*mEigVect)(i).assign_real (v (m));
277 
278  if (bPair)
279  {
280  (*mEigVect)(i).assign_imag (- v (m + 1));
281  m += 2;
282  bPair = false;
283  }
284  else
285  {
286  (*mEigVect)(i).assign_imag (v (m + 1));
287  bPair = true;
288  }
289  }
290  else
291  {
292  const float zero(0.F);
293  (*mEigVect)(i).assign_real (v (m));
294  (*mEigVect)(i).set_imag (zero);
295  m++;
296  }
297  }
298  }
299  else
300  {
301  lWork = -1;
302 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
303  SHSEQR (Chars::pE(), 1, Chars::pN(), 1,
304 #else
305  SHSEQR (Chars::pE(), Chars::pN(),
306 #endif
307  &nM, &ilo, &ihi, mA, &nM, vR, vI, nullptr, &nM, &dWork, &lWork, &nOutInfo);
308  _check_negative(CVM_WRONGMKLARG, nOutInfo);
309  lWork = static_cast<tint>(dWork) * 11;
310  if (lWork > vWork.size()) vWork.resize(lWork);
311 
312 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
313  SHSEQR (Chars::pE(), 1, Chars::pN(), 1,
314 #else
315  SHSEQR (Chars::pE(), Chars::pN(),
316 #endif
317  &nM, &ilo, &ihi, mA, &nM, vR, vI, nullptr, &nM, vWork, &lWork, &nOutInfo);
318  _check_negative(CVM_WRONGMKLARG, nOutInfo);
319  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SHSEQR", __FILE__, __LINE__);
320  }
321 
322  vRes.assign_real (vR);
323  vRes.assign_imag (vI);
324  }
325 }
326 
327 template<>
328 CVM_API void
329 __eig<basic_cvector<double, std::complex<double> >, basic_srmatrix<double>, basic_scmatrix<double, std::complex<double> > >
331  const basic_srmatrix<double>& mArg,
333  bool bRightVect) throw (cvmexception)
334 {
335  const bool bEigVect = (mEigVect != nullptr);
336  const tint nM = mArg.msize();
337  tint lWork = -1;
338  tint ilo = 0;
339  tint ihi = 0;
340  tint nOutInfo = 0;
341 
342  _check_ne(CVM_SIZESMISMATCH, vRes.size(), nM);
343  if (nM == 1)
344  {
345  vRes[CVM0] = std::complex<double>(mArg(CVM0,CVM0), 0.);
346  if (bEigVect)
347  {
348  const std::complex<double> one(1., 0.);
349  mEigVect -> resize (1);
350  (*mEigVect)[CVM0].set(one);
351  }
352  }
353  else
354  {
355  basic_srmatrix<double> mA (mArg);
356  basic_rvector<double> vScale (nM);
357  basic_rvector<double> vTau (_cvm_max<tint>(1, nM - 1));
358  basic_rvector<double> vR (nM);
359  basic_rvector<double> vI (nM);
360 
361  DGEBAL (Chars::pB(),
362 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
363  1,
364 #endif
365  &nM, mA, &nM, &ilo, &ihi, vScale, &nOutInfo);
366  _check_negative(CVM_WRONGMKLARG, nOutInfo);
367 
368  double dWork;
369  DGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, &dWork, &lWork, &nOutInfo);
370  _check_negative(CVM_WRONGMKLARG, nOutInfo);
371  lWork = static_cast<tint>(dWork);
372  basic_rvector<double> vWork (lWork);
373 
374  DGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, vWork, &lWork, &nOutInfo);
375  _check_negative(CVM_WRONGMKLARG, nOutInfo);
376 
377  if (bEigVect)
378  {
379  const tint one(1);
380  tint m = 0;
381  tint lSelect = 0;
382  const tint ldvl = bRightVect ? 1 : nM;
383  const tint ldvr = bRightVect ? nM : 1;
384  basic_srmatrix<double> vl (ldvl);
385  basic_srmatrix<double> vr (ldvr);
386  basic_rvector <double> work (3 * nM);
387  basic_srmatrix<double> mH (mA);
388  const char* pRL = bRightVect ? Chars::pR() : Chars::pL();
389 
390  if (bRightVect)
391  {
392  vr = mA;
393  }
394  else
395  {
396  vl = mA;
397  }
398 
399  lWork = -1;
400  DORGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, &dWork, &lWork, &nOutInfo);
401  _check_negative(CVM_WRONGMKLARG, nOutInfo);
402  lWork = static_cast<tint>(dWork);
403  if (lWork > vWork.size()) vWork.resize(lWork);
404 
405  DORGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, vWork, &lWork, &nOutInfo);
406  _check_negative(CVM_WRONGMKLARG, nOutInfo);
407 
408  lWork = -1;
409 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
410  DHSEQR (Chars::pS(), 1, Chars::pV(), 1,
411 #else
412  DHSEQR (Chars::pS(), Chars::pV(),
413 #endif
414  &nM, &ilo, &ihi, mH, &nM, vR, vI, bRightVect ? vr : vl, &nM, &dWork, &lWork, &nOutInfo);
415  _check_negative(CVM_WRONGMKLARG, nOutInfo);
416  lWork = static_cast<tint>(dWork) * 11;
417  if (lWork > vWork.size()) vWork.resize(lWork);
418 
419 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
420  DHSEQR (Chars::pS(), 1, Chars::pV(), 1,
421 #else
422  DHSEQR (Chars::pS(), Chars::pV(),
423 #endif
424  &nM, &ilo, &ihi, mH, &nM, vR, vI, bRightVect ? vr : vl, &nM, vWork, &lWork, &nOutInfo);
425  _check_negative(CVM_WRONGMKLARG, nOutInfo);
426  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DHSEQR", __FILE__, __LINE__);
427 
428 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
429  DTREVC (pRL, 1, Chars::pB(), 1,
430 #else
431  DTREVC (pRL, Chars::pB(),
432 #endif
433  &lSelect, &nM, mH, &nM, vl, &ldvl, vr, &ldvr, &nM, &m, work, &nOutInfo);
434  _check_negative(CVM_WRONGMKLARG, nOutInfo);
435 
436  basic_srmatrix<double>& v = bRightVect ? vr : vl;
437  const tint ldv = bRightVect ? ldvr : ldvl;
438 
439 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
440  DGEBAK (Chars::pB(), 1, pRL, 1,
441 #else
442  DGEBAK (Chars::pB(), pRL,
443 #endif
444  &nM, &ilo, &ihi, vScale, &nM, v, &ldv, &nOutInfo);
445  _check_negative(CVM_WRONGMKLARG, nOutInfo);
446 
447  bool bPair = false;
448  m = CVM0;
449  mEigVect -> resize (nM);
450  for (tint i = CVM0; i < nM + CVM0; i++)
451  {
452  if (_abs(vI(i)) > basic_cvmMachMin<double>())
453  {
454  (*mEigVect)(i).assign_real (v (m));
455 
456  if (bPair)
457  {
458  (*mEigVect)(i).assign_imag (- v (m + 1));
459  m += 2;
460  bPair = false;
461  }
462  else
463  {
464  (*mEigVect)(i).assign_imag (v (m + 1));
465  bPair = true;
466  }
467  }
468  else
469  {
470  const double zero(0.);
471  (*mEigVect)(i).assign_real (v (m));
472  (*mEigVect)(i).set_imag (zero);
473  m++;
474  }
475  }
476  }
477  else
478  {
479  lWork = -1;
480 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
481  DHSEQR (Chars::pE(), 1, Chars::pN(), 1,
482 #else
483  DHSEQR (Chars::pE(), Chars::pN(),
484 #endif
485  &nM, &ilo, &ihi, mA, &nM, vR, vI, nullptr, &nM, &dWork, &lWork, &nOutInfo);
486  _check_negative(CVM_WRONGMKLARG, nOutInfo);
487  lWork = static_cast<tint>(dWork) * 11;
488  if (lWork > vWork.size()) vWork.resize(lWork);
489 
490 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
491  DHSEQR (Chars::pE(), 1, Chars::pN(), 1,
492 #else
493  DHSEQR (Chars::pE(), Chars::pN(),
494 #endif
495  &nM, &ilo, &ihi, mA, &nM, vR, vI, nullptr, &nM, vWork, &lWork, &nOutInfo);
496  _check_negative(CVM_WRONGMKLARG, nOutInfo);
497  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DHSEQR", __FILE__, __LINE__);
498  }
499 
500  vRes.assign_real (vR);
501  vRes.assign_imag (vI);
502  }
503 }
504 
505 template<>
506 CVM_API void
507 __eig<basic_cvector<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> > >
509  const basic_scmatrix<float, std::complex<float> >& mArg,
510  basic_scmatrix<float, std::complex<float> >* mEigVect,
511  bool bRightVect) throw (cvmexception)
512 {
513  const bool bEigVect = (mEigVect != nullptr);
514  const tint nM = mArg.msize();
515  tint lWork = -1;
516  tint ilo = 0;
517  tint ihi = 0;
518  tint nOutInfo = 0;
519 
520  _check_ne(CVM_SIZESMISMATCH, vRes.size(), nM);
521  if (nM == 1)
522  {
523  vRes[CVM0] = mArg(CVM0,CVM0);
524  if (bEigVect)
525  {
526  const std::complex<float> one(1.F, 0.F);
527  mEigVect -> resize (1);
528  (*mEigVect)[CVM0].set(one);
529  }
530  }
531  else
532  {
533  basic_scmatrix<float, std::complex<float> > mA (mArg);
534  basic_rvector<float> vScale (nM);
535  basic_cvector<float, std::complex<float> > vTau (_cvm_max<tint>(1, nM - 1));
537 
538  CGEBAL (Chars::pB(),
539 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
540  1,
541 #endif
542  &nM, mA, &nM, &ilo, &ihi, vScale, &nOutInfo);
543  _check_negative(CVM_WRONGMKLARG, nOutInfo);
544 
545  std::complex<float> dWork;
546  CGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, &dWork, &lWork, &nOutInfo);
547  lWork = static_cast<tint>(dWork.real());
549 
550  CGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, vWork, &lWork, &nOutInfo);
551  _check_negative(CVM_WRONGMKLARG, nOutInfo);
552 
553  if (bEigVect)
554  {
555  const tint one(1);
556  tint m = 0;
557  tint lSelect = 0;
558  const tint ldvl = bRightVect ? 1 : nM;
559  const tint ldvr = bRightVect ? nM : 1;
560  basic_scmatrix<float, std::complex<float> > vl (ldvl);
561  basic_scmatrix<float, std::complex<float> > vr (ldvr);
563  basic_rvector <float> rwork (nM);
564  const char* pRL = bRightVect ? Chars::pR() : Chars::pL();
565 
566  if (bRightVect)
567  {
568  vr = mA;
569  }
570  else
571  {
572  vl = mA;
573  }
574 
575  lWork = -1;
576  CUNGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, &dWork, &lWork, &nOutInfo);
577  lWork = static_cast<tint>(dWork.real());
578  if (lWork > vWork.size()) vWork.resize(lWork);
579 
580  CUNGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, vWork, &lWork, &nOutInfo);
581  _check_negative(CVM_WRONGMKLARG, nOutInfo);
582 
583  lWork = -1;
584 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
585  CHSEQR (Chars::pS(), 1, Chars::pV(), 1,
586 #else
587  CHSEQR (Chars::pS(), Chars::pV(),
588 #endif
589  &nM, &ilo, &ihi, mA, &nM, vW, bRightVect ? vr : vl, &nM, &dWork, &lWork, &nOutInfo);
590  lWork = static_cast<tint> (dWork.real()) * 11;
591  if (lWork > vWork.size()) vWork.resize(lWork);
592 
593 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
594  CHSEQR (Chars::pS(), 1, Chars::pV(), 1,
595 #else
596  CHSEQR (Chars::pS(), Chars::pV(),
597 #endif
598  &nM, &ilo, &ihi, mA, &nM, vW, bRightVect ? vr : vl, &nM, vWork, &lWork, &nOutInfo);
599  _check_negative(CVM_WRONGMKLARG, nOutInfo);
600  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CHSEQR", __FILE__, __LINE__);
601 
602 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
603  CTREVC (pRL, 1, Chars::pB(), 1,
604 #else
605  CTREVC (pRL, Chars::pB(),
606 #endif
607  &lSelect, &nM, mA, &nM, vl, &ldvl, vr, &ldvr, &nM, &m, work, rwork, &nOutInfo);
608  _check_negative(CVM_WRONGMKLARG, nOutInfo);
609 
610  basic_scmatrix<float, std::complex<float> >& v = bRightVect ? vr : vl;
611  const tint ldv = bRightVect ? ldvr : ldvl;
612 
613 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
614  CGEBAK (Chars::pB(), 1, pRL, 1,
615 #else
616  CGEBAK (Chars::pB(), pRL,
617 #endif
618  &nM, &ilo, &ihi, vScale, &nM, v, &ldv, &nOutInfo);
619  _check_negative(CVM_WRONGMKLARG, nOutInfo);
620 
621  (*mEigVect) << v;
622  }
623  else
624  {
625  lWork = -1;
626 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
627  CHSEQR (Chars::pE(), 1, Chars::pN(), 1,
628 #else
629  CHSEQR (Chars::pE(), Chars::pN(),
630 #endif
631  &nM, &ilo, &ihi, mA, &nM, vW, nullptr, &nM, &dWork, &lWork, &nOutInfo);
632  lWork = static_cast<tint> (dWork.real()) * 11;
633  if (lWork > vWork.size()) vWork.resize(lWork);
634 
635 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
636  CHSEQR (Chars::pE(), 1, Chars::pN(), 1,
637 #else
638  CHSEQR (Chars::pE(), Chars::pN(),
639 #endif
640  &nM, &ilo, &ihi, mA, &nM, vW, nullptr, &nM, vWork, &lWork, &nOutInfo);
641  _check_negative(CVM_WRONGMKLARG, nOutInfo);
642  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CHSEQR", __FILE__, __LINE__);
643  }
644 
645  vRes.assign (vW, vW.incr());
646  }
647 }
648 
649 template<>
650 CVM_API void
651 __eig<basic_cvector<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> > >
653  const basic_scmatrix<double, std::complex<double> >& mArg,
654  basic_scmatrix<double, std::complex<double> >* mEigVect,
655  bool bRightVect) throw (cvmexception)
656 {
657  const bool bEigVect = (mEigVect != nullptr);
658  const tint nM = mArg.msize();
659  tint lWork = -1;
660  tint ilo = 0;
661  tint ihi = 0;
662  tint nOutInfo = 0;
663 
664  _check_ne(CVM_SIZESMISMATCH, vRes.size(), nM);
665  if (nM == 1)
666  {
667  vRes[CVM0] = mArg(CVM0,CVM0);
668  if (bEigVect)
669  {
670  const std::complex<double> one(1., 0.);
671  mEigVect -> resize (1);
672  (*mEigVect)[CVM0].set(one);
673  }
674  }
675  else
676  {
677  basic_scmatrix<double, std::complex<double> > mA (mArg);
678  basic_rvector<double> vScale (nM);
679  basic_cvector<double, std::complex<double> > vTau (_cvm_max<tint>(1, nM - 1));
681 
682  ZGEBAL (Chars::pB(),
683 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
684  1,
685 #endif
686  &nM, mA, &nM, &ilo, &ihi, vScale, &nOutInfo);
687  _check_negative(CVM_WRONGMKLARG, nOutInfo);
688 
689  std::complex<double> dWork;
690  ZGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, &dWork, &lWork, &nOutInfo);
691  lWork = static_cast<tint> (dWork.real());
693 
694  ZGEHRD (&nM, &ilo, &ihi, mA, &nM, vTau, vWork, &lWork, &nOutInfo);
695  _check_negative(CVM_WRONGMKLARG, nOutInfo);
696 
697  if (bEigVect)
698  {
699  const tint one(1);
700  tint m = 0;
701  tint lSelect = 0;
702  const tint ldvl = bRightVect ? 1 : nM;
703  const tint ldvr = bRightVect ? nM : 1;
704  basic_scmatrix<double, std::complex<double> > vl (ldvl);
705  basic_scmatrix<double, std::complex<double> > vr (ldvr);
707  basic_rvector <double> rwork (nM);
708  const char* pRL = bRightVect ? Chars::pR() : Chars::pL();
709 
710  if (bRightVect)
711  {
712  vr = mA;
713  }
714  else
715  {
716  vl = mA;
717  }
718 
719  lWork = -1;
720  ZUNGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, &dWork, &lWork, &nOutInfo);
721  lWork = static_cast<tint> (dWork.real());
722  if (lWork > vWork.size()) vWork.resize(lWork);
723 
724  ZUNGHR (&nM, &one, &nM, bRightVect ? vr : vl, &nM, vTau, vWork, &lWork, &nOutInfo);
725  _check_negative(CVM_WRONGMKLARG, nOutInfo);
726 
727  lWork = -1;
728 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
729  ZHSEQR (Chars::pS(), 1, Chars::pV(), 1,
730 #else
731  ZHSEQR (Chars::pS(), Chars::pV(),
732 #endif
733  &nM, &ilo, &ihi, mA, &nM, vW, bRightVect ? vr : vl, &nM, &dWork, &lWork, &nOutInfo);
734  lWork = static_cast<tint> (dWork.real()) * 11;
735  if (lWork > vWork.size()) vWork.resize(lWork);
736 
737 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
738  ZHSEQR (Chars::pS(), 1, Chars::pV(), 1,
739 #else
740  ZHSEQR (Chars::pS(), Chars::pV(),
741 #endif
742  &nM, &ilo, &ihi, mA, &nM, vW, bRightVect ? vr : vl, &nM, vWork, &lWork, &nOutInfo);
743  _check_negative(CVM_WRONGMKLARG, nOutInfo);
744  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZHSEQR", __FILE__, __LINE__);
745 
746 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
747  ZTREVC (pRL, 1, Chars::pB(), 1,
748 #else
749  ZTREVC (pRL, Chars::pB(),
750 #endif
751  &lSelect, &nM, mA, &nM, vl, &ldvl, vr, &ldvr, &nM, &m, work, rwork, &nOutInfo);
752  _check_negative(CVM_WRONGMKLARG, nOutInfo);
753 
754  basic_scmatrix<double, std::complex<double> >& v = bRightVect ? vr : vl;
755  const tint ldv = bRightVect ? ldvr : ldvl;
756 
757 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
758  ZGEBAK (Chars::pB(), 1, pRL, 1,
759 #else
760  ZGEBAK (Chars::pB(), pRL,
761 #endif
762  &nM, &ilo, &ihi, vScale, &nM, v, &ldv, &nOutInfo);
763  _check_negative(CVM_WRONGMKLARG, nOutInfo);
764 
765  (*mEigVect) << v;
766  }
767  else
768  {
769  lWork = -1;
770 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
771  ZHSEQR (Chars::pE(), 1, Chars::pN(), 1,
772 #else
773  ZHSEQR (Chars::pE(), Chars::pN(),
774 #endif
775  &nM, &ilo, &ihi, mA, &nM, vW, nullptr, &nM, &dWork, &lWork, &nOutInfo);
776  lWork = static_cast<tint> (dWork.real()) * 11;
777  if (lWork > vWork.size()) vWork.resize(lWork);
778 
779 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
780  ZHSEQR (Chars::pE(), 1, Chars::pN(), 1,
781 #else
782  ZHSEQR (Chars::pE(), Chars::pN(),
783 #endif
784  &nM, &ilo, &ihi, mA, &nM, vW, nullptr, &nM, vWork, &lWork, &nOutInfo);
785  _check_negative(CVM_WRONGMKLARG, nOutInfo);
786  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZHSEQR", __FILE__, __LINE__);
787  }
788 
789  vRes.assign (vW, vW.incr());
790  }
791 }
792 
793 // internal helper
794 void _harvest_eig_vect(basic_scmatrix<float, std::complex<float> >& mEigVect,
795  const basic_srmatrix<float>& v,
796  const basic_rvector<float>& alphai)
797 {
798  static float zero(0.F);
799  for (tint i = CVM0; i < v.nsize() + CVM0; ++i) {
800  mEigVect(i).assign_real(v(i));
801  if (i < v.nsize() - (TINT_ONE - CVM0) && alphai(i) > basic_cvmMachMin<treal>()) { // (i, i+1) form conjugate pair
802  mEigVect(i + 1).assign_real(v(i));
803  mEigVect(i).assign_imag(v(i + 1));
804  mEigVect(i + 1).assign_imag(-v(i + 1));
805  ++i;
806  continue;
807  }
808  mEigVect(i).set_imag(zero);
809  }
810 }
811 
812 void _harvest_eig_vect(basic_scmatrix<double, std::complex<double> >& mEigVect,
813  const basic_srmatrix<double>& v,
814  const basic_rvector<double>& alphai)
815 {
816  static double zero(0.);
817  for (tint i = CVM0; i < v.nsize() + CVM0; ++i) {
818  mEigVect(i).assign_real(v(i));
819  if (i < v.nsize() - (TINT_ONE - CVM0) && alphai(i) > basic_cvmMachMin<treal>()) { // (i, i+1) form conjugate pair
820  mEigVect(i + 1).assign_real(v(i));
821  mEigVect(i).assign_imag(v(i + 1));
822  mEigVect(i + 1).assign_imag(-v(i + 1));
823  ++i;
824  continue;
825  }
826  mEigVect(i).set_imag(zero);
827  }
828 }
829 
830 // Computes the generalized eigenvalues, and the
831 // left and / or right generalized eigenvectors for a pair
832 // of nonsymmetric matrices.
833 // Overrides A and B!
834 template<>
835 CVM_API void
836 __ggev<basic_srmatrix<float>, basic_scmatrix<float, std::complex<float> >,
838 (basic_srmatrix<float>& mA,
839 basic_srmatrix<float>& mB,
841 basic_rvector<float>& vBeta,
842 basic_scmatrix<float, std::complex<float> >* mEigVectLeft,
843 basic_scmatrix<float, std::complex<float> >* mEigVectRight) throw (cvmexception)
844 {
845  const char* jobvl = mEigVectLeft == nullptr ? Chars::pN() : Chars::pV();
846  const char* jobvr = mEigVectRight == nullptr ? Chars::pN() : Chars::pV();
847  const tint n = mA.nsize(); // assert equal sizes outside
848  basic_rvector<float> alphar(n);
849  basic_rvector<float> alphai(n);
850  basic_srmatrix<float> vl(n);
851  basic_srmatrix<float> vr(n);
852  tint lWork = -1;
853  float dWork;
854  tint nOutInfo = 0;
855 
856  SGGEV(jobvl,
857 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
858  1,
859 #endif
860  jobvr,
861 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
862  1,
863 #endif
864  &n, mA, mA._pld(), mB, mB._pld(),
865  alphar, alphai, vBeta,
866  vl, vl._pld(), vr, vr._pld(),
867  &dWork, &lWork, &nOutInfo);
868  _check_negative(CVM_WRONGMKLARG, nOutInfo);
869 
870  lWork = static_cast<tint>(dWork);
871  basic_rvector<float> vWork(lWork);
872  SGGEV(jobvl,
873 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
874  1,
875 #endif
876  jobvr,
877 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
878  1,
879 #endif
880  &n, mA, mA._pld(), mB, mB._pld(),
881  alphar, alphai, vBeta,
882  vl, vl._pld(), vr, vr._pld(),
883  vWork, vWork._psize(), &nOutInfo);
884  _check_negative(CVM_WRONGMKLARG, nOutInfo);
885 
886  if (mEigVectLeft != nullptr) {
887  _harvest_eig_vect(*mEigVectLeft, vl, alphai);
888  }
889  if (mEigVectRight != nullptr) {
890  _harvest_eig_vect(*mEigVectRight, vr, alphai);
891  }
892  vAlpha = basic_cvector<float, std::complex<float> >(alphar, alphai);
893 }
894 
895 
896 template<>
897 CVM_API void
898 __ggev<basic_srmatrix<double>, basic_scmatrix<double, std::complex<double> >,
900  (basic_srmatrix<double>& mA,
901  basic_srmatrix<double>& mB,
903  basic_rvector<double>& vBeta,
904  basic_scmatrix<double, std::complex<double> >* mEigVectLeft,
905  basic_scmatrix<double, std::complex<double> >* mEigVectRight) throw (cvmexception)
906 {
907  const char* jobvl = mEigVectLeft == nullptr ? Chars::pN() : Chars::pV();
908  const char* jobvr = mEigVectRight == nullptr ? Chars::pN() : Chars::pV();
909  const tint n = mA.nsize(); // assert equal sizes outside
910  basic_rvector<double> alphar(n);
911  basic_rvector<double> alphai(n);
912  basic_srmatrix<double> vl(n);
913  basic_srmatrix<double> vr(n);
914  tint lWork = -1;
915  double dWork;
916  tint nOutInfo = 0;
917 
918  DGGEV(jobvl,
919 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
920  1,
921 #endif
922  jobvr,
923 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
924  1,
925 #endif
926  &n, mA, mA._pld(), mB, mB._pld(),
927  alphar, alphai, vBeta,
928  vl, vl._pld(), vr, vr._pld(),
929  &dWork, &lWork, &nOutInfo);
930  _check_negative(CVM_WRONGMKLARG, nOutInfo);
931 
932  lWork = static_cast<tint>(dWork);
933  basic_rvector<double> vWork(lWork);
934  DGGEV(jobvl,
935 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
936  1,
937 #endif
938  jobvr,
939 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
940  1,
941 #endif
942  &n, mA, mA._pld(), mB, mB._pld(),
943  alphar, alphai, vBeta,
944  vl, vl._pld(), vr, vr._pld(),
945  vWork, vWork._psize(), &nOutInfo);
946  _check_negative(CVM_WRONGMKLARG, nOutInfo);
947 
948  if (mEigVectLeft != nullptr) {
949  _harvest_eig_vect(*mEigVectLeft, vl, alphai);
950  }
951  if (mEigVectRight != nullptr) {
952  _harvest_eig_vect(*mEigVectRight, vr, alphai);
953  }
954  vAlpha = basic_cvector<double, std::complex<double> >(alphar, alphai);
955 }
956 
957 template<>
958 CVM_API void
959 __ggev<basic_scmatrix<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> >,
960  basic_cvector<float, std::complex<float> >, basic_cvector<float, std::complex<float> > >
961  (basic_scmatrix<float, std::complex<float> >& mA,
962  basic_scmatrix<float, std::complex<float> >& mB,
963  basic_cvector<float, std::complex<float> >& vAlpha,
964  basic_cvector<float, std::complex<float> >& vBeta,
965  basic_scmatrix<float, std::complex<float> >* mEigVectLeft,
966  basic_scmatrix<float, std::complex<float> >* mEigVectRight) throw (cvmexception)
967 {
968  const char* jobvl = mEigVectLeft == nullptr ? Chars::pN() : Chars::pV();
969  const char* jobvr = mEigVectRight == nullptr ? Chars::pN() : Chars::pV();
970  const tint n = mA.nsize(); // assert equal sizes outside
971  basic_scmatrix<float, std::complex<float> > vl(n);
972  basic_scmatrix<float, std::complex<float> > vr(n);
973  basic_rvector<float> vRWork(8 * n);
974  tint lWork = -1;
975  std::complex<float> dWork;
976  tint nOutInfo = 0;
977 
978  CGGEV(jobvl,
979 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
980  1,
981 #endif
982  jobvr,
983 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
984  1,
985 #endif
986  &n, mA, mA._pld(), mB, mB._pld(),
987  vAlpha, vBeta,
988  vl, vl._pld(), vr, vr._pld(),
989  &dWork, &lWork, vRWork, &nOutInfo);
990  _check_negative(CVM_WRONGMKLARG, nOutInfo);
991 
992  lWork = static_cast<tint>(dWork.real());
993  basic_cvector<float, std::complex<float> > vWork(lWork);
994  CGGEV(jobvl,
995 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
996  1,
997 #endif
998  jobvr,
999 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1000  1,
1001 #endif
1002  &n, mA, mA._pld(), mB, mB._pld(),
1003  vAlpha, vBeta,
1004  vl, vl._pld(), vr, vr._pld(),
1005  vWork, vWork._psize(), vRWork, &nOutInfo);
1006  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1007 
1008  if (mEigVectLeft != nullptr) {
1009  *mEigVectLeft = vl;
1010  }
1011  if (mEigVectRight != nullptr) {
1012  *mEigVectRight = vr;
1013  }
1014 }
1015 
1016 template<>
1017 CVM_API void
1018 __ggev<basic_scmatrix<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> >,
1019  basic_cvector<double, std::complex<double> >, basic_cvector<double, std::complex<double> > >
1020  (basic_scmatrix<double, std::complex<double> >& mA,
1021  basic_scmatrix<double, std::complex<double> >& mB,
1022  basic_cvector<double, std::complex<double> >& vAlpha,
1023  basic_cvector<double, std::complex<double> >& vBeta,
1024  basic_scmatrix<double, std::complex<double> >* mEigVectLeft,
1025  basic_scmatrix<double, std::complex<double> >* mEigVectRight) throw (cvmexception)
1026 {
1027  const char* jobvl = mEigVectLeft == nullptr ? Chars::pN() : Chars::pV();
1028  const char* jobvr = mEigVectRight == nullptr ? Chars::pN() : Chars::pV();
1029  const tint n = mA.nsize(); // assert equal sizes outside
1030  basic_scmatrix<double, std::complex<double> > vl(n);
1031  basic_scmatrix<double, std::complex<double> > vr(n);
1032  basic_rvector<double> vRWork(8 * n);
1033  tint lWork = -1;
1034  std::complex<double> dWork;
1035  tint nOutInfo = 0;
1036 
1037  ZGGEV(jobvl,
1038 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1039  1,
1040 #endif
1041  jobvr,
1042 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1043  1,
1044 #endif
1045  &n, mA, mA._pld(), mB, mB._pld(),
1046  vAlpha, vBeta,
1047  vl, vl._pld(), vr, vr._pld(),
1048  &dWork, &lWork, vRWork, &nOutInfo);
1049  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1050 
1051  lWork = static_cast<tint>(dWork.real());
1052  basic_cvector<double, std::complex<double> > vWork(lWork);
1053  ZGGEV(jobvl,
1054 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1055  1,
1056 #endif
1057  jobvr,
1058 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1059  1,
1060 #endif
1061  &n, mA, mA._pld(), mB, mB._pld(),
1062  vAlpha, vBeta,
1063  vl, vl._pld(), vr, vr._pld(),
1064  vWork, vWork._psize(), vRWork, &nOutInfo);
1065  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1066 
1067  if (mEigVectLeft != nullptr) {
1068  *mEigVectLeft = vl;
1069  }
1070  if (mEigVectRight != nullptr) {
1071  *mEigVectRight = vr;
1072  }
1073 }
1074 
1076