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
rvector.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 
12 CVM_NAMESPACE_BEG
13 
14 template<>
15 CVM_API float __dot<float> (const float* mpd, tint mnSize, tint mnIncr, const float* pd, tint nIncr)
16 {
17  return SDOT (&mnSize, mpd, &mnIncr, pd, &nIncr);
18 }
19 
20 template<>
21 CVM_API double __dot<double> (const double* mpd, tint mnSize, tint mnIncr, const double* pd, tint nIncr)
22 {
23  return DDOT (&mnSize, mpd, &mnIncr, pd, &nIncr);
24 }
25 
26 template<>
27 CVM_API void
28 __gemv<float, basic_rmatrix<float>, basic_rvector<float> >
29  (bool bLeft,
30  const basic_rmatrix<float>& m,
31  float dAlpha,
32  const basic_rvector<float>& v,
33  float dBeta,
34  basic_rvector<float>& vRes)
35 {
36  SGEMV (bLeft ? Chars::pT() : Chars::pN(),
37 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
38  1,
39 #endif
40  m._pm(), m._pn(), &dAlpha, m._pd(), m._pldm(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
41 }
42 
43 template<>
44 CVM_API void
45 __gemv<double, basic_rmatrix<double>, basic_rvector<double> >
46  (bool bLeft,
47  const basic_rmatrix<double>& m,
48  double dAlpha,
49  const basic_rvector<double>& v,
50  double dBeta,
51  basic_rvector<double>& vRes)
52 {
53  DGEMV (bLeft ? Chars::pT() : Chars::pN(),
54 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
55  1,
56 #endif
57  m._pm(), m._pn(), &dAlpha, m._pd(), m._pldm(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
58 }
59 
60 template<>
61 CVM_API void
62 __gbmv<float, basic_srbmatrix<float>, basic_rvector<float> >
63  (bool bLeft,
64  const basic_srbmatrix<float>& m,
65  float dAlpha,
66  const basic_rvector<float>& v,
67  float dBeta,
68  basic_rvector<float>& vRes)
69 {
70  SGBMV (bLeft ? Chars::pT() : Chars::pN(),
71 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
72  1,
73 #endif
74  m._pm(), m._pn(), m._pl(), m._pu(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
75 }
76 
77 template<>
78 CVM_API void
79 __gbmv<double, basic_srbmatrix<double>, basic_rvector<double> >
80  (bool bLeft,
81  const basic_srbmatrix<double>& m,
82  double dAlpha,
83  const basic_rvector<double>& v,
84  double dBeta,
85  basic_rvector<double>& vRes)
86 {
87  DGBMV (bLeft ? Chars::pT() : Chars::pN(),
88 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
89  1,
90 #endif
91  m._pm(), m._pn(), m._pl(), m._pu(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
92 }
93 
94 template<>
95 CVM_API void
96 __symv<float, basic_srsmatrix<float>, basic_rvector<float> >
97  (const basic_srsmatrix<float>& m,
98  float dAlpha,
99  const basic_rvector<float>& v,
100  float dBeta,
101  basic_rvector<float>& vRes)
102 {
103  SSYMV (Chars::pU(),
104 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
105  1,
106 #endif
107  m._pm(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
108 }
109 
110 template<>
111 CVM_API void
112 __symv<double, basic_srsmatrix<double>, basic_rvector<double> >
113  (const basic_srsmatrix<double>& m,
114  double dAlpha,
115  const basic_rvector<double>& v,
116  double dBeta,
117  basic_rvector<double>& vRes)
118 {
119  DSYMV (Chars::pU(),
120 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
121  1,
122 #endif
123  m._pm(), &dAlpha, m, m._pld(), v, v._pincr(), &dBeta, vRes, vRes._pincr());
124 }
125 
126 template<>
127 CVM_API void
128 __svd<float, basic_rmatrix<float>, basic_srmatrix<float> >
129  (float* pd, tint nSize, tint nIncr,
130  const basic_rmatrix<float>& mArg,
131  basic_srmatrix<float>* mU,
132  basic_srmatrix<float>* mVH) throw (cvmexception)
133 {
134  const bool bSimple = (mU == nullptr || mVH == nullptr);
135  const tint nM = mArg.msize();
136  const tint nN = mArg.nsize();
137  const tint m = _cvm_min<tint>(nM, nN);
138  const tint m4 = m * 4;
139  const tint zero(0);
140  const tint one(1);
141  tint lwork = -1; // to calculate lwork
142  tint nOutInfo = 0;
143 
144  _check_ne(CVM_SIZESMISMATCH, m, nSize);
145 
146  basic_rvector<float> mD (nSize);
147  mD.assign (pd, nIncr);
148 
149  basic_rvector<float> vOffDiag (_cvm_max<tint>(1, m - 1));
150  basic_rvector<float> vTauQ (m);
151  basic_rvector<float> vTauP (m);
152  basic_rmatrix<float> mA (mArg);
153  float dwork;
154 
155  SGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
156  lwork = static_cast<tint>(dwork);
157  basic_rvector<float> vwork (lwork);
158  SGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
159  _check_negative(CVM_WRONGMKLARG, nOutInfo);
160 
161  if (bSimple)
162  {
163  if (m4 > vwork.size()) vwork.resize(m4);
164  SBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
165 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
166  1,
167 #endif
168  &m,
169  &zero, &zero, &zero,
170  mD, vOffDiag,
171  nullptr, &one, nullptr, &one, nullptr, &one,
172  vwork, &nOutInfo);
173  }
174  else
175  {
176  basic_rmatrix<float> Q (mA);
177  basic_rmatrix<float> P (mA);
178 
179  if (nM > nN) Q.resize(nM, nM);
180  if (nM < nN) P.resize(nN, nN);
181 
182  lwork = -1;
183  SORGBR (Chars::pQ(),
184 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
185  1,
186 #endif
187  &nM, &nM, &nN,
188  Q, Q._pld(),
189  vTauQ,
190  &dwork, &lwork, &nOutInfo);
191  _check_negative(CVM_WRONGMKLARG, nOutInfo);
192  lwork = static_cast<tint>(dwork);
193  if (lwork > vwork.size()) vwork.resize(lwork);
194 
195  SORGBR (Chars::pQ(),
196 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
197  1,
198 #endif
199  &nM, &nM, &nN,
200  Q, Q._pld(),
201  vTauQ,
202  vwork, &lwork, &nOutInfo);
203  _check_negative(CVM_WRONGMKLARG, nOutInfo);
204 
205  lwork = -1;
206  SORGBR (Chars::pP(),
207 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
208  1,
209 #endif
210  &nN, &nN, &nM,
211  P, P._pld(),
212  vTauP,
213  &dwork, &lwork, &nOutInfo);
214  _check_negative(CVM_WRONGMKLARG, nOutInfo);
215  lwork = static_cast<tint>(dwork);
216  if (lwork > vwork.size()) vwork.resize(lwork);
217 
218  SORGBR (Chars::pP(),
219 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
220  1,
221 #endif
222  &nN, &nN, &nM,
223  P, P._pld(),
224  vTauP,
225  vwork, &lwork, &nOutInfo);
226  _check_negative(CVM_WRONGMKLARG, nOutInfo);
227 
228  if (m4 > vwork.size()) vwork.resize(m4);
229  SBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
230 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
231  1,
232 #endif
233  &m,
234  &nN, &nM,
235  &zero,
236  mD, vOffDiag,
237  P, P._pld(), Q, Q._pld(), nullptr, &one,
238  vwork, &nOutInfo);
239 
240  (*mU) = basic_srmatrix<float>(Q.resize(nM, nM));
241  (*mVH) = basic_srmatrix<float>(P.resize(nN, nN));
242  }
243  _check_negative(CVM_WRONGMKLARG, nOutInfo);
244  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SBDSQR", __FILE__, __LINE__);
245 
246  __copy<float> (nSize, mD, mD.incr(), pd, nIncr);
247 }
248 
249 template<>
250 CVM_API void
251 __svd<double, basic_rmatrix<double>, basic_srmatrix<double> >
252  (double* pd, tint nSize, tint nIncr,
253  const basic_rmatrix<double>& mArg,
254  basic_srmatrix<double>* mU,
255  basic_srmatrix<double>* mVH) throw (cvmexception)
256 {
257  const bool bSimple = (mU == nullptr || mVH == nullptr);
258  const tint nM = mArg.msize();
259  const tint nN = mArg.nsize();
260  const tint m = _cvm_min<tint>(nM, nN);
261  const tint m4 = m * 4;
262  const tint zero(0);
263  const tint one(1);
264  tint lwork = -1; // to calculate lwork
265  tint nOutInfo = 0;
266 
267  _check_ne(CVM_SIZESMISMATCH, m, nSize);
268 
269  basic_rvector<double> mD (nSize);
270  mD.assign (pd, nIncr);
271 
272  basic_rvector<double> vOffDiag (_cvm_max<tint>(1, m - 1));
273  basic_rvector<double> vTauQ (m);
274  basic_rvector<double> vTauP (m);
275  basic_rmatrix<double> mA (mArg);
276  double dwork;
277 
278  DGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
279  lwork = static_cast<tint>(dwork);
280  basic_rvector<double> vwork(lwork);
281  DGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
282  _check_negative(CVM_WRONGMKLARG, nOutInfo);
283 
284  if (bSimple)
285  {
286  if (m4 > vwork.size()) vwork.resize(m4);
287  DBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
288 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
289  1,
290 #endif
291  &m,
292  &zero, &zero, &zero,
293  mD, vOffDiag,
294  nullptr, &one, nullptr, &one, nullptr, &one,
295  vwork, &nOutInfo);
296  }
297  else
298  {
299  basic_rmatrix<double> Q (mA);
300  basic_rmatrix<double> P (mA);
301 
302  if (nM > nN) Q.resize(nM, nM);
303  if (nM < nN) P.resize(nN, nN);
304 
305  lwork = -1;
306  DORGBR (Chars::pQ(),
307 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
308  1,
309 #endif
310  &nM, &nM, &nN,
311  Q, Q._pld(),
312  vTauQ,
313  &dwork, &lwork, &nOutInfo);
314  _check_negative(CVM_WRONGMKLARG, nOutInfo);
315  lwork = static_cast<tint>(dwork);
316  if (lwork > vwork.size()) vwork.resize(lwork);
317 
318  DORGBR (Chars::pQ(),
319 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
320  1,
321 #endif
322  &nM, &nM, &nN,
323  Q, Q._pld(),
324  vTauQ,
325  vwork, &lwork, &nOutInfo);
326  _check_negative(CVM_WRONGMKLARG, nOutInfo);
327 
328  lwork = -1;
329  DORGBR (Chars::pP(),
330 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
331  1,
332 #endif
333  &nN, &nN, &nM,
334  P, P._pld(),
335  vTauP,
336  &dwork, &lwork, &nOutInfo);
337  _check_negative(CVM_WRONGMKLARG, nOutInfo);
338  lwork = static_cast<tint>(dwork);
339  if (lwork > vwork.size()) vwork.resize(lwork);
340 
341  DORGBR (Chars::pP(),
342 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
343  1,
344 #endif
345  &nN, &nN, &nM,
346  P, P._pld(),
347  vTauP,
348  vwork, &lwork, &nOutInfo);
349  _check_negative(CVM_WRONGMKLARG, nOutInfo);
350 
351  if (m4 > vwork.size()) vwork.resize(m4);
352  DBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
353 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
354  1,
355 #endif
356  &m,
357  &nN, &nM,
358  &zero,
359  mD, vOffDiag,
360  P, P._pld(), Q, Q._pld(), nullptr, &one,
361  vwork, &nOutInfo);
362 
363  (*mU) = basic_srmatrix<double>(Q.resize(nM, nM));
364  (*mVH) = basic_srmatrix<double>(P.resize(nN, nN));
365  }
366  _check_negative(CVM_WRONGMKLARG, nOutInfo);
367  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DBDSQR", __FILE__, __LINE__);
368 
369  __copy<double> (nSize, mD, mD.incr(), pd, nIncr);
370 }
371 
373 template<>
374 CVM_API void
375 __svd<float, basic_cmatrix<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> > >
376  (float* pd, tint nSize, tint nIncr,
377  const basic_cmatrix<float, std::complex<float> >& mArg,
378  basic_scmatrix<float, std::complex<float> >* mU,
379  basic_scmatrix<float, std::complex<float> >* mVH) throw (cvmexception)
380 {
381  const bool bSimple = (mU == nullptr || mVH == nullptr);
382  const tint nM = mArg.msize();
383  const tint nN = mArg.nsize();
384  const tint m = _cvm_min<tint>(nM, nN);
385  const tint zero(0);
386  const tint one(1);
387  tint lwork = -1; // to calculate lwork
388  tint nOutInfo = 0;
389 
390  _check_ne(CVM_SIZESMISMATCH, m, nSize);
391 
392  basic_rvector<float> mD (nSize);
393  mD.assign (pd, nIncr);
394 
395  basic_rvector<float> vOffDiag (_cvm_max<tint>(1, m - 1));
396  basic_cvector<float, std::complex<float> > vTauQ (m);
397  basic_cvector<float, std::complex<float> > vTauP (m);
398  basic_cmatrix<float, std::complex<float> > mA (mArg);
399  basic_rvector<float> rwork (m * 4);
400  std::complex<float> dwork;
401 
402  CGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
403  lwork = static_cast<tint>(dwork.real());
404  basic_cvector<float, std::complex<float> > vwork (lwork);
405  CGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
406  _check_negative(CVM_WRONGMKLARG, nOutInfo);
407 
408  if (bSimple)
409  {
410  CBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
411 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
412  1,
413 #endif
414  &m,
415  &zero, &zero, &zero,
416  mD, vOffDiag,
417  nullptr, &one, nullptr, &one, nullptr, &one,
418  rwork, &nOutInfo);
419  }
420  else
421  {
422  basic_cmatrix<float, std::complex<float> > Q (mA);
423  basic_cmatrix<float, std::complex<float> > P (mA);
424 
425  if (nM > nN) Q.resize(nM, nM);
426  if (nM < nN) P.resize(nN, nN);
427 
428  lwork = -1;
429  CUNGBR (Chars::pQ(),
430 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
431  1,
432 #endif
433  &nM, &nM, &nN,
434  Q, Q._pld(),
435  vTauQ,
436  &dwork, &lwork, &nOutInfo);
437  _check_negative(CVM_WRONGMKLARG, nOutInfo);
438  lwork = static_cast<tint>(dwork.real());
439  if (lwork > vwork.size()) vwork.resize(lwork);
440 
441  CUNGBR (Chars::pQ(),
442 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
443  1,
444 #endif
445  &nM, &nM, &nN,
446  Q, Q._pld(),
447  vTauQ,
448  vwork, &lwork, &nOutInfo);
449  _check_negative(CVM_WRONGMKLARG, nOutInfo);
450 
451  lwork = -1;
452  CUNGBR (Chars::pP(),
453 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
454  1,
455 #endif
456  &nN, &nN, &nM,
457  P, P._pld(),
458  vTauP,
459  &dwork, &lwork, &nOutInfo);
460  _check_negative(CVM_WRONGMKLARG, nOutInfo);
461  lwork = static_cast<tint>(dwork.real());
462  if (lwork > vwork.size()) vwork.resize(lwork);
463 
464  CUNGBR (Chars::pP(),
465 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
466  1,
467 #endif
468  &nN, &nN, &nM,
469  P, P._pld(),
470  vTauP,
471  vwork, &lwork, &nOutInfo);
472  _check_negative(CVM_WRONGMKLARG, nOutInfo);
473 
474  CBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
475 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
476  1,
477 #endif
478  &m,
479  &nN, &nM,
480  &zero,
481  mD, vOffDiag,
482  P, P._pld(), Q, Q._pld(), nullptr, &one,
483  rwork, &nOutInfo);
484 
485  (*mU) = basic_scmatrix<float, std::complex<float> >(Q.resize(nM, nM));
486  (*mVH) = basic_scmatrix<float, std::complex<float> >(P.resize(nN, nN));
487  }
488  _check_negative(CVM_WRONGMKLARG, nOutInfo);
489  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CBDSQR", __FILE__, __LINE__);
490 
491  __copy<float> (nSize, mD, mD.incr(), pd, nIncr);
492 }
493 
494 template<>
495 CVM_API void
496 __svd<double, basic_cmatrix<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> > >
497  (double* pd, tint nSize, tint nIncr,
498  const basic_cmatrix<double, std::complex<double> >& mArg,
499  basic_scmatrix<double, std::complex<double> >* mU,
500  basic_scmatrix<double, std::complex<double> >* mVH) throw (cvmexception)
501 {
502  const bool bSimple = (mU == nullptr || mVH == nullptr);
503  const tint nM = mArg.msize();
504  const tint nN = mArg.nsize();
505  const tint m = _cvm_min<tint>(nM, nN);
506  const tint zero(0);
507  const tint one(1);
508  tint lwork = -1; // to calculate lwork
509  tint nOutInfo = 0;
510 
511  _check_ne(CVM_SIZESMISMATCH, m, nSize);
512 
513  basic_rvector<double> mD (nSize);
514  mD.assign (pd, nIncr);
515 
516  basic_rvector<double> vOffDiag (_cvm_max<tint>(1, m - 1));
517  basic_cvector<double, std::complex<double> > vTauQ (m);
518  basic_cvector<double, std::complex<double> > vTauP (m);
519  basic_cmatrix<double, std::complex<double> > mA (mArg);
520  basic_rvector<double> rwork (m * 4);
521  std::complex<double> dwork;
522 
523  ZGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
524  lwork = static_cast<tint>(dwork.real());
525  basic_cvector<double, std::complex<double> > vwork (lwork);
526  ZGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
527  _check_negative(CVM_WRONGMKLARG, nOutInfo);
528 
529  if (bSimple)
530  {
531  ZBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
532 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
533  1,
534 #endif
535  &m,
536  &zero, &zero, &zero,
537  mD, vOffDiag,
538  nullptr, &one, nullptr, &one, nullptr, &one,
539  rwork, &nOutInfo);
540  }
541  else
542  {
543  basic_cmatrix<double, std::complex<double> > Q (mA);
544  basic_cmatrix<double, std::complex<double> > P (mA);
545 
546  if (nM > nN) Q.resize(nM, nM);
547  if (nM < nN) P.resize(nN, nN);
548 
549  lwork = -1;
550  ZUNGBR (Chars::pQ(),
551 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
552  1,
553 #endif
554  &nM, &nM, &nN,
555  Q, Q._pld(),
556  vTauQ,
557  &dwork, &lwork, &nOutInfo);
558  _check_negative(CVM_WRONGMKLARG, nOutInfo);
559  lwork = static_cast<tint>(dwork.real());
560  if (lwork > vwork.size()) vwork.resize(lwork);
561 
562  ZUNGBR (Chars::pQ(),
563 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
564  1,
565 #endif
566  &nM, &nM, &nN,
567  Q, Q._pld(),
568  vTauQ,
569  vwork, &lwork, &nOutInfo);
570  _check_negative(CVM_WRONGMKLARG, nOutInfo);
571 
572  lwork = -1;
573  ZUNGBR (Chars::pP(),
574 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
575  1,
576 #endif
577  &nN, &nN, &nM,
578  P, P._pld(),
579  vTauP,
580  &dwork, &lwork, &nOutInfo);
581  _check_negative(CVM_WRONGMKLARG, nOutInfo);
582  lwork = static_cast<tint>(dwork.real());
583  if (lwork > vwork.size()) vwork.resize(lwork);
584 
585  ZUNGBR (Chars::pP(),
586 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
587  1,
588 #endif
589  &nN, &nN, &nM,
590  P, P._pld(),
591  vTauP,
592  vwork, &lwork, &nOutInfo);
593  _check_negative(CVM_WRONGMKLARG, nOutInfo);
594 
595  ZBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
596 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
597  1,
598 #endif
599  &m,
600  &nN, &nM,
601  &zero,
602  mD, vOffDiag,
603  P, P._pld(), Q, Q._pld(), nullptr, &one,
604  rwork, &nOutInfo);
605 
606  (*mU) = basic_scmatrix<double, std::complex<double> >(Q.resize(nM, nM));
607  (*mVH) = basic_scmatrix<double, std::complex<double> >(P.resize(nN, nN));
608  }
609  _check_negative(CVM_WRONGMKLARG, nOutInfo);
610  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZBDSQR", __FILE__, __LINE__);
611 
612  __copy<double> (nSize, mD, mD.incr(), pd, nIncr);
613 }
615 
616 template<>
617 CVM_API void
618 __svd<float, basic_srbmatrix<float>, basic_srmatrix<float> >
619  (float* pd, tint nSize, tint nIncr,
620  const basic_srbmatrix<float>& mArg,
621  basic_srmatrix<float>* mU,
622  basic_srmatrix<float>* mVH) throw (cvmexception)
623 {
624  const tint zero(0);
625  const tint m = mArg.msize();
626 
627  _check_ne(CVM_SIZESMISMATCH, m, nSize);
628 
629  const bool bSimple = (mU == nullptr || mVH == nullptr);
630  tint nOutInfo = 0;
631 
632  basic_rvector<float> mD (nSize);
633  mD.assign (pd, nIncr);
634 
635  basic_rvector <float> vOffDiag (_cvm_max<tint>(1, m - 1));
636  basic_srmatrix <float> mQ (bSimple ? 1 : m);
637  basic_srmatrix <float> mPT (bSimple ? 1 : m);
638  basic_srmatrix <float> mC (1);
639  basic_rvector <float> vwork (m * 4); // 2*max(M,N), but M==N, then sharing it with 4*N
640  basic_srbmatrix<float> mA (mArg);
641 
642  SGBBRD (bSimple ? Chars::pN() : Chars::pB(),
643 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
644  1,
645 #endif
646  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
647  mD, vOffDiag,
648  mQ, mQ._pm(),
649  mPT, mPT._pm(),
650  mC, mC._pm(),
651  vwork, &nOutInfo);
652  _check_negative(CVM_WRONGMKLARG, nOutInfo);
653 
654  SBDSQR (Chars::pU(),
655 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
656  1,
657 #endif
658  &m,
659  bSimple ? &zero : &m,
660  bSimple ? &zero : &m,
661  &zero,
662  mD, vOffDiag,
663  mPT, mPT._pm(),
664  mQ, mQ._pm(),
665  mC, mC._pm(),
666  vwork, &nOutInfo);
667  _check_negative(CVM_WRONGMKLARG, nOutInfo);
668  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SBDSQR", __FILE__, __LINE__);
669 
670  if (!bSimple)
671  {
672  (*mU) = mQ;
673  (*mVH) = mPT;
674  }
675 
676  __copy<float> (nSize, mD, mD.incr(), pd, nIncr);
677 }
678 
679 template<>
680 CVM_API void
681 __svd<double, basic_srbmatrix<double>, basic_srmatrix<double> >
682  (double* pd, tint nSize, tint nIncr,
683  const basic_srbmatrix<double>& mArg,
684  basic_srmatrix<double>* mU,
685  basic_srmatrix<double>* mVH) throw (cvmexception)
686 {
687  const tint zero(0);
688  const tint m = mArg.msize();
689 
690  _check_ne(CVM_SIZESMISMATCH, m, nSize);
691 
692  const bool bSimple = (mU == nullptr || mVH == nullptr);
693  tint nOutInfo = 0;
694 
695  basic_rvector<double> mD (nSize);
696  mD.assign (pd, nIncr);
697 
698  basic_rvector <double> vOffDiag (_cvm_max<tint>(1, m - 1));
699  basic_srmatrix <double> mQ (bSimple ? 1 : m);
700  basic_srmatrix <double> mPT (bSimple ? 1 : m);
701  basic_srmatrix <double> mC (1);
702  basic_rvector <double> vwork (m * 4); // 2*max(M,N), but M==N, so sharing it with 4*N
703  basic_srbmatrix<double> mA (mArg);
704 
705  DGBBRD (bSimple ? Chars::pN() : Chars::pB(),
706 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
707  1,
708 #endif
709  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
710  mD, vOffDiag,
711  mQ, mQ._pm(),
712  mPT, mPT._pm(),
713  mC, mC._pm(),
714  vwork, &nOutInfo);
715  _check_negative(CVM_WRONGMKLARG, nOutInfo);
716 
717  DBDSQR (Chars::pU(),
718 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
719  1,
720 #endif
721  &m,
722  bSimple ? &zero : &m,
723  bSimple ? &zero : &m,
724  &zero,
725  mD, vOffDiag,
726  mPT, mPT._pm(),
727  mQ, mQ._pm(),
728  mC, mC._pm(),
729  vwork, &nOutInfo);
730  _check_negative(CVM_WRONGMKLARG, nOutInfo);
731  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DBDSQR", __FILE__, __LINE__);
732 
733  if (!bSimple)
734  {
735  (*mU) = mQ;
736  (*mVH) = mPT;
737  }
738 
739  __copy<double> (nSize, mD, mD.incr(), pd, nIncr);
740 }
741 
743 template<>
744 CVM_API void
745 __svd<float, basic_scbmatrix<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> > >
746  (float* pd, tint nSize, tint nIncr,
747  const basic_scbmatrix<float, std::complex<float> >& mArg,
748  basic_scmatrix<float, std::complex<float> >* mU,
749  basic_scmatrix<float, std::complex<float> >* mVH) throw (cvmexception)
750 {
751  const tint zero(0);
752  const tint m = mArg.msize();
753 
754  _check_ne(CVM_SIZESMISMATCH, m, nSize);
755 
756  const bool bSimple = (mU == nullptr || mVH == nullptr);
757  tint nOutInfo = 0;
758 
759  basic_rvector<float> mD (nSize);
760  mD.assign (pd, nIncr);
761 
762  basic_rvector <float> vOffDiag (_cvm_max<tint>(1, m - 1));
763  basic_scmatrix <float, std::complex<float> > mQ (bSimple ? 1 : m);
764  basic_scmatrix <float, std::complex<float> > mPT (bSimple ? 1 : m);
765  basic_scmatrix <float, std::complex<float> > mC (1);
766  basic_cvector <float, std::complex<float> > vwork (m);
767  basic_scbmatrix<float, std::complex<float> > mA (mArg);
768  basic_rvector <float> vrwork (m * 4); // max(M,N), but M==N, so sharing it with 4*N
769 
770  CGBBRD (bSimple ? Chars::pN() : Chars::pB(),
771 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
772  1,
773 #endif
774  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
775  mD, vOffDiag,
776  mQ, mQ._pm(),
777  mPT, mPT._pm(),
778  mC, mC._pm(),
779  vwork, vrwork, &nOutInfo);
780  _check_negative(CVM_WRONGMKLARG, nOutInfo);
781 
782  CBDSQR (Chars::pU(),
783 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
784  1,
785 #endif
786  &m,
787  bSimple ? &zero : &m,
788  bSimple ? &zero : &m,
789  &zero,
790  mD, vOffDiag,
791  mPT, mPT._pm(),
792  mQ, mQ._pm(),
793  mC, mC._pm(),
794  vrwork, &nOutInfo);
795  _check_negative(CVM_WRONGMKLARG, nOutInfo);
796  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CBDSQR", __FILE__, __LINE__);
797 
798  if (!bSimple)
799  {
800  (*mU) = mQ;
801  (*mVH) = mPT;
802  }
803 
804  __copy<float> (nSize, mD, mD.incr(), pd, nIncr);
805 }
806 
807 template<>
808 CVM_API void
809 __svd<double, basic_scbmatrix<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> > >
810  (double* pd, tint nSize, tint nIncr,
811  const basic_scbmatrix<double, std::complex<double> >& mArg,
812  basic_scmatrix<double, std::complex<double> >* mU,
813  basic_scmatrix<double, std::complex<double> >* mVH) throw (cvmexception)
814 {
815  const tint zero(0);
816  const tint m = mArg.msize();
817 
818  _check_ne(CVM_SIZESMISMATCH, m, nSize);
819 
820  const bool bSimple = (mU == nullptr || mVH == nullptr);
821  tint nOutInfo = 0;
822 
823  basic_rvector<double> mD (nSize);
824  mD.assign (pd, nIncr);
825 
826  basic_rvector <double> vOffDiag (_cvm_max<tint>(1, m - 1));
827  basic_scmatrix <double, std::complex<double> > mQ (bSimple ? 1 : m);
828  basic_scmatrix <double, std::complex<double> > mPT (bSimple ? 1 : m);
829  basic_scmatrix <double, std::complex<double> > mC (1);
830  basic_cvector <double, std::complex<double> > vwork (m);
831  basic_scbmatrix<double, std::complex<double> > mA (mArg);
832  basic_rvector <double> vrwork (m * 4); // max(M,N), but M==N, so sharing it with 4*N
833 
834  ZGBBRD (bSimple ? Chars::pN() : Chars::pB(),
835 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
836  1,
837 #endif
838  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
839  mD, vOffDiag,
840  mQ, mQ._pm(),
841  mPT, mPT._pm(),
842  mC, mC._pm(),
843  vwork, vrwork, &nOutInfo);
844  _check_negative(CVM_WRONGMKLARG, nOutInfo);
845 
846  basic_rvector<double> vwork2 (m * 4);
847  ZBDSQR (Chars::pU(),
848 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
849  1,
850 #endif
851  &m,
852  bSimple ? &zero : &m,
853  bSimple ? &zero : &m,
854  &zero,
855  mD, vOffDiag,
856  mPT, mPT._pm(),
857  mQ, mQ._pm(),
858  mC, mC._pm(),
859  vrwork, &nOutInfo);
860  _check_negative(CVM_WRONGMKLARG, nOutInfo);
861  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZBDSQR", __FILE__, __LINE__);
862 
863  if (!bSimple)
864  {
865  (*mU) = mQ;
866  (*mVH) = mPT;
867  }
868 
869  __copy<double> (nSize, mD, mD.incr(), pd, nIncr);
870 }
872 
873 template<>
874 CVM_API void
875 __eig<basic_rvector<float>, basic_srsmatrix<float>, basic_srmatrix<float> >
876  (basic_rvector<float>& vRes,
877  const basic_srsmatrix<float>& mArg,
878  basic_srmatrix<float>* mEigVect,
879  bool /*bRightVect*/) throw (cvmexception)
880 {
881  const tint nM = mArg.msize();
882  _check_ne(CVM_SIZESMISMATCH, nM, vRes.size());
883  const bool bEigVect = (mEigVect != nullptr);
884 
885  if (nM == 1)
886  {
887  vRes[CVM0] = mArg(CVM0,CVM0);
888  if (bEigVect)
889  {
890  const float one(1.F);
891  mEigVect -> resize (1);
892  (*mEigVect)[CVM0].set (one);
893  }
894  }
895  else
896  {
897  const char* pcJob = bEigVect ? Chars::pV() : Chars::pN();
898  basic_srsmatrix<float> mA (mArg);
899  tint lwork = -1, liwork = -1;
900  float work_size = 0.;
901  tint iwork_size = 0;
902  tint nOutInfo = 0;
903 
904  SSYEVD (pcJob,
905  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
906  1,
907  #endif
908  Chars::pU(),
909  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
910  1,
911  #endif
912  &nM, mA, mA._pld(), vRes, &work_size, &lwork, &iwork_size, &liwork, &nOutInfo);
913  lwork = static_cast<tint> (work_size);
914  liwork = iwork_size;
915  basic_rvector<float> work (lwork);
916  basic_array<tint,tint> iwork (liwork);
917 
918  SSYEVD (pcJob,
919  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
920  1,
921  #endif
922  Chars::pU(),
923  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
924  1,
925  #endif
926  &nM, mA, mA._pld(), vRes, work, &lwork, iwork, &liwork, &nOutInfo);
927  _check_negative(CVM_WRONGMKLARG, nOutInfo);
928  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SSYEVD", __FILE__, __LINE__);
929 
930  if (bEigVect)
931  {
932  (*mEigVect) << mA;
933  }
934  }
935 }
936 
937 template<>
938 CVM_API void
939 __eig<basic_rvector<double>, basic_srsmatrix<double>, basic_srmatrix<double> >
940  (basic_rvector<double>& vRes,
941  const basic_srsmatrix<double>& mArg,
942  basic_srmatrix<double>* mEigVect,
943  bool /*bRightVect*/) throw (cvmexception)
944 {
945  const tint nM = mArg.msize();
946  _check_ne(CVM_SIZESMISMATCH, nM, vRes.size());
947  const bool bEigVect = (mEigVect != nullptr);
948 
949  if (nM == 1)
950  {
951  vRes[CVM0] = mArg(CVM0,CVM0);
952  if (bEigVect)
953  {
954  const double one(1.);
955  mEigVect -> resize (1);
956  (*mEigVect)[CVM0].set (one);
957  }
958  }
959  else
960  {
961  const char* pcJob = bEigVect ? Chars::pV() : Chars::pN();
962  basic_srsmatrix<double> mA (mArg);
963  tint lwork = -1, liwork = -1;
964  double work_size = 0.;
965  tint iwork_size = 0;
966  tint nOutInfo = 0;
967 
968  DSYEVD (pcJob,
969  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
970  1,
971  #endif
972  Chars::pU(),
973  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
974  1,
975  #endif
976  &nM, mA, mA._pld(), vRes, &work_size, &lwork, &iwork_size, &liwork, &nOutInfo);
977  lwork = static_cast<tint> (work_size);
978  liwork = iwork_size;
979  basic_rvector<double> work (lwork);
980  basic_array<tint,tint> iwork (liwork);
981 
982  DSYEVD (pcJob,
983  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
984  1,
985  #endif
986  Chars::pU(),
987  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
988  1,
989  #endif
990  &nM, mA, mA._pld(), vRes, work, &lwork, iwork, &liwork, &nOutInfo);
991  _check_negative(CVM_WRONGMKLARG, nOutInfo);
992  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DSYEVD", __FILE__, __LINE__);
993 
994  if (bEigVect)
995  {
996  (*mEigVect) << mA;
997  }
998  }
999 }
1000 
1002 template<>
1003 CVM_API void
1004 __eig<basic_rvector<float>, basic_schmatrix<float, std::complex<float> >, basic_scmatrix<float, std::complex<float> > >
1005  (basic_rvector<float>& vRes,
1006  const basic_schmatrix<float, std::complex<float> >& mArg,
1007  basic_scmatrix<float, std::complex<float> >* mEigVect,
1008  bool /*bRightVect*/) throw (cvmexception)
1009 {
1010  const tint nM = mArg.msize();
1011  _check_ne(CVM_SIZESMISMATCH, nM, vRes.size());
1012  const bool bEigVect = (mEigVect != nullptr);
1013 
1014  if (nM == 1)
1015  {
1016  vRes[CVM0] = 1.F;
1017  if (bEigVect)
1018  {
1019  mEigVect -> resize (1);
1020  (*mEigVect)[CVM0].set (mArg(CVM0,CVM0));
1021  }
1022  }
1023  else
1024  {
1025  const char* pcJob = bEigVect ? Chars::pV() : Chars::pN();
1026  basic_schmatrix<float, std::complex<float> > mA (mArg);
1027  tint lwork = -1, lrwork = -1, liwork = -1;
1028  std::complex<float> work_size;
1029  float rwork_size = 0.;
1030  tint iwork_size = 0;
1031  tint nOutInfo = 0;
1032 
1033  CHEEVD (pcJob,
1034  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1035  1,
1036  #endif
1037  Chars::pU(),
1038  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1039  1,
1040  #endif
1041  &nM, mA, mA._pld(), vRes, &work_size, &lwork, &rwork_size, &lrwork, &iwork_size, &liwork, &nOutInfo);
1042  lwork = static_cast<tint> (work_size.real());
1043  lrwork = static_cast<tint> (rwork_size);
1044  liwork = iwork_size;
1045  basic_cvector<float, std::complex<float> > work (lwork);
1046  basic_rvector<float> rwork (lrwork);
1047  basic_array<tint,tint> iwork (liwork);
1048 
1049  CHEEVD (pcJob,
1050  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1051  1,
1052  #endif
1053  Chars::pU(),
1054  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1055  1,
1056  #endif
1057  &nM, mA, mA._pld(), vRes, work, &lwork, rwork, &lrwork, iwork, &liwork, &nOutInfo);
1058  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1059  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CHEEVD", __FILE__, __LINE__);
1060 
1061  if (bEigVect)
1062  {
1063  (*mEigVect) << mA;
1064  }
1065  }
1066 }
1067 
1068 template<>
1069 CVM_API void
1070 __eig<basic_rvector<double>, basic_schmatrix<double, std::complex<double> >, basic_scmatrix<double, std::complex<double> > >
1071  (basic_rvector<double>& vRes,
1072  const basic_schmatrix<double, std::complex<double> >& mArg,
1073  basic_scmatrix<double, std::complex<double> >* mEigVect,
1074  bool /*bRightVect*/) throw (cvmexception)
1075 {
1076  const tint nM = mArg.msize();
1077  _check_ne(CVM_SIZESMISMATCH, nM, vRes.size());
1078  const bool bEigVect = (mEigVect != nullptr);
1079 
1080  if (nM == 1)
1081  {
1082  vRes[CVM0] = 1.;
1083  if (bEigVect)
1084  {
1085  mEigVect -> resize (1);
1086  (*mEigVect)[CVM0].set (mArg(CVM0,CVM0));
1087  }
1088  }
1089  else
1090  {
1091  const char* pcJob = bEigVect ? Chars::pV() : Chars::pN();
1092  basic_schmatrix<double, std::complex<double> > mA (mArg);
1093  tint lwork = -1, lrwork = -1, liwork = -1;
1094  std::complex<double> work_size;
1095  double rwork_size = 0.;
1096  tint iwork_size = 0;
1097  tint nOutInfo = 0;
1098 
1099  ZHEEVD (pcJob,
1100  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1101  1,
1102  #endif
1103  Chars::pU(),
1104  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1105  1,
1106  #endif
1107  &nM, mA, mA._pld(), vRes, &work_size, &lwork, &rwork_size, &lrwork, &iwork_size, &liwork, &nOutInfo);
1108  lwork = static_cast<tint> (work_size.real());
1109  lrwork = static_cast<tint> (rwork_size);
1110  liwork = iwork_size;
1111  basic_cvector<double, std::complex<double> > work (lwork);
1112  basic_rvector<double> rwork (lrwork);
1113  basic_array<tint,tint> iwork (liwork);
1114 
1115  ZHEEVD (pcJob,
1116  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1117  1,
1118  #endif
1119  Chars::pU(),
1120  #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1121  1,
1122  #endif
1123  &nM, mA, mA._pld(), vRes, work, &lwork, rwork, &lrwork, iwork, &liwork, &nOutInfo);
1124  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1125  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZHEEVD", __FILE__, __LINE__);
1126 
1127  if (bEigVect)
1128  {
1129  (*mEigVect) << mA;
1130  }
1131  }
1132 }
1134 
1135 // pseudo (generalized) inversion routines
1136 template<>
1137 CVM_API void
1138 __pinv<float, basic_rmatrix<float>, basic_rmatrix<float> >
1139  (basic_rmatrix<float>& mX,
1140  const basic_rmatrix<float>& mArg, float threshold) throw (cvmexception)
1141 {
1142  const tint nM = mArg.msize();
1143  const tint nN = mArg.nsize();
1144  const tint m = _cvm_min<tint>(nM, nN);
1145  tint lwork = -1; // to calculate lwork
1146  tint nOutInfo = 0;
1147 
1148  basic_rvector<float> mD (m);
1149  basic_rvector<float> vOffDiag (_cvm_max<tint>(1, m - 1));
1150  basic_rvector<float> vTauQ (m);
1151  basic_rvector<float> vTauP (m);
1152  basic_rmatrix<float> mA (mArg);
1153  float dwork;
1154 
1155  SGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
1156  lwork = static_cast<tint>(dwork);
1157  basic_rvector<float> vwork(lwork);
1158  SGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
1159  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1160 
1161  const tint zero(0);
1162  const tint one(1);
1163 
1164  // few words about economy:
1165  // for m > n case we care about m-by-n matrix U and n-by-n matrix VH
1166  // for m < n case we care about m-by-m matrix U and m-by-n matrix VH
1167  // however, the whole matrix A is needed to start computations
1168  basic_rmatrix<float> Q (mA);
1169  basic_rmatrix<float> P (mA);
1170 
1171  {
1172  lwork = -1;
1173  SORGBR (Chars::pQ(),
1174 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1175  1,
1176 #endif
1177  &nM, &m, &nN,
1178  Q, Q._pld(),
1179  vTauQ,
1180  &dwork, &lwork, &nOutInfo);
1181  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1182 
1183  lwork = static_cast<tint>(dwork);
1184  basic_rvector<float> vwork3(lwork);
1185  SORGBR (Chars::pQ(),
1186 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1187  1,
1188 #endif
1189  &nM, &m, &nN,
1190  Q, Q._pld(),
1191  vTauQ,
1192  vwork3, &lwork, &nOutInfo);
1193  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1194  }
1195 
1196  {
1197  lwork = -1;
1198  SORGBR (Chars::pP(),
1199 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1200  1,
1201 #endif
1202  &m, &nN, &nM,
1203  P, P._pld(),
1204  vTauP,
1205  &dwork, &lwork, &nOutInfo);
1206  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1207 
1208  lwork = static_cast<tint>(dwork);
1209  basic_rvector<float> vwork3(lwork);
1210  SORGBR (Chars::pP(),
1211 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1212  1,
1213 #endif
1214  &m, &nN, &nM,
1215  P, P._pld(),
1216  vTauP,
1217  vwork3, &lwork, &nOutInfo);
1218  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1219  }
1220 
1221  basic_rvector<float> vwork2 (m * 4);
1222  SBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
1223 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1224  1,
1225 #endif
1226  &m,
1227  &nN, &nM,
1228  &zero,
1229  mD, vOffDiag,
1230  P, P._pld(), Q, Q._pld(), nullptr, &one,
1231  vwork2, &nOutInfo);
1232  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1233  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SBDSQR", __FILE__, __LINE__);
1234 
1235  if (nM > nN) P.resize(nN, nN); // VH
1236  if (nM < nN) Q.resize(nM, nM); // U
1237  for (tint i = CVM0; i < P.msize() + CVM0; ++i) {
1238  if (mD[i] > threshold) {
1239  P[i] /= mD[i]; // multiplying V by S^-1
1240  }
1241  else {
1242  P[i].set(0.F);
1243  }
1244  }
1245  mX.mult (~P, ~Q);
1246 }
1247 
1248 template<>
1249 CVM_API void
1250 __pinv<double, basic_rmatrix<double>, basic_rmatrix<double> >
1251  (basic_rmatrix<double>& mX,
1252  const basic_rmatrix<double>& mArg, double threshold) throw (cvmexception)
1253 {
1254  const tint nM = mArg.msize();
1255  const tint nN = mArg.nsize();
1256  const tint m = _cvm_min<tint>(nM, nN);
1257  tint lwork = -1; // to calculate lwork
1258  tint nOutInfo = 0;
1259 
1260  basic_rvector<double> mD (m);
1261  basic_rvector<double> vOffDiag (_cvm_max<tint>(1, m - 1));
1262  basic_rvector<double> vTauQ (m);
1263  basic_rvector<double> vTauP (m);
1264  basic_rmatrix<double> mA (mArg);
1265  double dwork;
1266 
1267  DGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
1268  lwork = static_cast<tint>(dwork);
1269  basic_rvector<double> vwork(lwork);
1270  DGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
1271  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1272 
1273  const tint zero(0);
1274  const tint one(1);
1275 
1276  // few words about economy:
1277  // for m > n case we care about m-by-n matrix U and n-by-n matrix VH
1278  // for m < n case we care about m-by-m matrix U and m-by-n matrix VH
1279  // however, the whole matrix A is needed to start computations
1280  basic_rmatrix<double> Q (mA);
1281  basic_rmatrix<double> P (mA);
1282 
1283  {
1284  lwork = -1;
1285  DORGBR (Chars::pQ(),
1286 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1287  1,
1288 #endif
1289  &nM, &m, &nN,
1290  Q, Q._pld(),
1291  vTauQ,
1292  &dwork, &lwork, &nOutInfo);
1293  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1294 
1295 
1296  lwork = static_cast<tint>(dwork);
1297  basic_rvector<double> vwork3(lwork);
1298  DORGBR (Chars::pQ(),
1299 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1300  1,
1301 #endif
1302  &nM, &m, &nN,
1303  Q, Q._pld(),
1304  vTauQ,
1305  vwork3, &lwork, &nOutInfo);
1306  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1307  }
1308 
1309  {
1310  lwork = -1;
1311  DORGBR (Chars::pP(),
1312 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1313  1,
1314 #endif
1315  &m, &nN, &nM,
1316  P, P._pld(),
1317  vTauP,
1318  &dwork, &lwork, &nOutInfo);
1319  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1320 
1321  lwork = static_cast<tint>(dwork);
1322  basic_rvector<double> vwork3(lwork);
1323  DORGBR (Chars::pP(),
1324 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1325  1,
1326 #endif
1327  &m, &nN, &nM,
1328  P, P._pld(),
1329  vTauP,
1330  vwork3, &lwork, &nOutInfo);
1331  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1332  }
1333 
1334  basic_rvector<double> vwork2 (m * 4);
1335  DBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
1336 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1337  1,
1338 #endif
1339  &m,
1340  &nN, &nM,
1341  &zero,
1342  mD, vOffDiag,
1343  P, P._pld(), Q, Q._pld(), nullptr, &one,
1344  vwork2, &nOutInfo);
1345  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1346  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DBDSQR", __FILE__, __LINE__);
1347 
1348  if (nM > nN) P.resize(nN, nN); // VH
1349  if (nM < nN) Q.resize(nM, nM); // U
1350  for (tint i = CVM0; i < P.msize() + CVM0; ++i) {
1351  if (mD[i] > threshold) {
1352  P[i] /= mD[i]; // multiplying V by S^-1
1353  }
1354  else {
1355  P[i].set(0.);
1356  }
1357  }
1358  mX.mult (~P, ~Q);
1359 }
1360 
1362 template<>
1363 CVM_API void
1364 __pinv<float, basic_cmatrix<float, std::complex<float> >, basic_cmatrix<float, std::complex<float> > >
1365  (basic_cmatrix<float, std::complex<float> >& mX,
1366  const basic_cmatrix<float, std::complex<float> >& mArg, float threshold) throw (cvmexception)
1367 {
1368  const tint nM = mArg.msize();
1369  const tint nN = mArg.nsize();
1370  const tint m = _cvm_min<tint>(nM, nN);
1371  tint lwork = -1; // to calculate lwork
1372  tint nOutInfo = 0;
1373 
1374  basic_rvector<float> mD (m);
1375  basic_rvector<float> vOffDiag (_cvm_max<tint>(1, m - 1));
1376  basic_cvector<float, std::complex<float> > vTauQ (m);
1377  basic_cvector<float, std::complex<float> > vTauP (m);
1378  basic_cmatrix<float, std::complex<float> > mA (mArg);
1379  std::complex<float> dwork;
1380 
1381  CGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
1382  lwork = static_cast<tint>(dwork.real());
1383  basic_cvector<float, std::complex<float> > vwork(lwork);
1384  CGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
1385  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1386 
1387  const tint zero(0);
1388  const tint one(1);
1389 
1390  // few words about economy:
1391  // for m > n case we care about m-by-n matrix U and n-by-n matrix VH
1392  // for m < n case we care about m-by-m matrix U and m-by-n matrix VH
1393  // however, the whole matrix A is needed to start computations
1394  basic_cmatrix<float, std::complex<float> > Q (mA);
1395  basic_cmatrix<float, std::complex<float> > P (mA);
1396 
1397  {
1398  lwork = -1;
1399  CUNGBR (Chars::pQ(),
1400 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1401  1,
1402 #endif
1403  &nM, &m, &nN,
1404  Q, Q._pld(),
1405  vTauQ,
1406  &dwork, &lwork, &nOutInfo);
1407  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1408 
1409  lwork = static_cast<tint>(dwork.real());
1410  basic_cvector<float, std::complex<float> > vwork3 (lwork);
1411  CUNGBR (Chars::pQ(),
1412 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1413  1,
1414 #endif
1415  &nM, &m, &nN,
1416  Q, Q._pld(),
1417  vTauQ,
1418  vwork3, &lwork, &nOutInfo);
1419  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1420  }
1421 
1422  {
1423  lwork = -1;
1424  CUNGBR (Chars::pP(),
1425 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1426  1,
1427 #endif
1428  &m, &nN, &nM,
1429  P, P._pld(),
1430  vTauP,
1431  &dwork, &lwork, &nOutInfo);
1432  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1433 
1434  lwork = static_cast<tint>(dwork.real());
1435  basic_cvector<float, std::complex<float> > vwork3 (lwork);
1436  CUNGBR (Chars::pP(),
1437 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1438  1,
1439 #endif
1440  &m, &nN, &nM,
1441  P, P._pld(),
1442  vTauP,
1443  vwork3, &lwork, &nOutInfo);
1444  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1445  }
1446 
1447  basic_rvector<float> vwork2 (m * 4);
1448  CBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
1449 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1450  1,
1451 #endif
1452  &m,
1453  &nN, &nM,
1454  &zero,
1455  mD, vOffDiag,
1456  P, P._pld(), Q, Q._pld(), nullptr, &one,
1457  vwork2, &nOutInfo);
1458  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1459  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CBDSQR", __FILE__, __LINE__);
1460 
1461  if (nM > nN) P.resize(nN, nN); // VH
1462  if (nM < nN) Q.resize(nM, nM); // U
1463  for (tint i = CVM0; i < P.msize() + CVM0; ++i) {
1464  if (mD[i] > threshold) {
1465  P[i] /= mD[i]; // multiplying V by S^-1
1466  }
1467  else {
1468  P[i].set(0.F);
1469  }
1470  }
1471  mX.mult (~P, ~Q);
1472 }
1473 
1474 template<>
1475 CVM_API void
1476 __pinv<double, basic_cmatrix<double, std::complex<double> >, basic_cmatrix<double, std::complex<double> > >
1477  (basic_cmatrix<double, std::complex<double> >& mX,
1478  const basic_cmatrix<double, std::complex<double> >& mArg, double threshold) throw (cvmexception)
1479 {
1480  const tint nM = mArg.msize();
1481  const tint nN = mArg.nsize();
1482  const tint m = _cvm_min<tint>(nM, nN);
1483  tint lwork = -1; // to calculate lwork
1484  tint nOutInfo = 0;
1485 
1486  basic_rvector<double> mD (m);
1487  basic_rvector<double> vOffDiag (_cvm_max<tint>(1, m - 1));
1488  basic_cvector<double, std::complex<double> > vTauQ (m);
1489  basic_cvector<double, std::complex<double> > vTauP (m);
1490  basic_cmatrix<double, std::complex<double> > mA (mArg);
1491  std::complex<double> dwork;
1492 
1493  ZGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, &dwork, &lwork, &nOutInfo);
1494  lwork = static_cast<tint>(dwork.real());
1495  basic_cvector<double, std::complex<double> > vwork(lwork);
1496  ZGEBRD (&nM, &nN, mA, mA._pld(), mD, vOffDiag, vTauQ, vTauP, vwork, &lwork, &nOutInfo);
1497  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1498 
1499  const tint zero(0);
1500  const tint one(1);
1501 
1502  // few words about economy:
1503  // for m > n case we care about m-by-n matrix U and n-by-n matrix VH
1504  // for m < n case we care about m-by-m matrix U and m-by-n matrix VH
1505  // however, the whole matrix A is needed to start computations
1506  basic_cmatrix<double, std::complex<double> > Q (mA);
1507  basic_cmatrix<double, std::complex<double> > P (mA);
1508 
1509  {
1510  lwork = -1;
1511  ZUNGBR (Chars::pQ(),
1512 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1513  1,
1514 #endif
1515  &nM, &m, &nN,
1516  Q, Q._pld(),
1517  vTauQ,
1518  &dwork, &lwork, &nOutInfo);
1519  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1520 
1521  lwork = static_cast<tint>(dwork.real());
1522  basic_cvector<double, std::complex<double> > vwork3 (lwork);
1523  ZUNGBR (Chars::pQ(),
1524 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1525  1,
1526 #endif
1527  &nM, &m, &nN,
1528  Q, Q._pld(),
1529  vTauQ,
1530  vwork3, &lwork, &nOutInfo);
1531  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1532  }
1533 
1534  {
1535  lwork = -1;
1536  ZUNGBR (Chars::pP(),
1537 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1538  1,
1539 #endif
1540  &m, &nN, &nM,
1541  P, P._pld(),
1542  vTauP,
1543  &dwork, &lwork, &nOutInfo);
1544  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1545 
1546  lwork = static_cast<tint>(dwork.real());
1547  basic_cvector<double, std::complex<double> > vwork3 (lwork);
1548  ZUNGBR (Chars::pP(),
1549 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1550  1,
1551 #endif
1552  &m, &nN, &nM,
1553  P, P._pld(),
1554  vTauP,
1555  vwork3, &lwork, &nOutInfo);
1556  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1557  }
1558 
1559  basic_rvector<double> vwork2 (m * 4);
1560  ZBDSQR (nM >= nN ? Chars::pU() : Chars::pL(),
1561 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1562  1,
1563 #endif
1564  &m,
1565  &nN, &nM,
1566  &zero,
1567  mD, vOffDiag,
1568  P, P._pld(), Q, Q._pld(), nullptr, &one,
1569  vwork2, &nOutInfo);
1570  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1571  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZBDSQR", __FILE__, __LINE__);
1572 
1573  if (nM > nN) P.resize(nN, nN); // VH
1574  if (nM < nN) Q.resize(nM, nM); // U
1575  for (tint i = CVM0; i < P.msize() + CVM0; ++i) {
1576  if (mD[i] > threshold) {
1577  P[i] /= mD[i]; // multiplying V by S^-1
1578  }
1579  else {
1580  P[i].set(0.);
1581  }
1582  }
1583  mX.mult (~P, ~Q);
1584 }
1586 
1587 template<>
1588 CVM_API void
1589 __pinv<float, basic_srbmatrix<float>, basic_rmatrix<float> >
1590  (basic_rmatrix<float>& mX,
1591  const basic_srbmatrix<float>& mArg, float threshold) throw (cvmexception)
1592 {
1593  const tint zero(0);
1594  const tint m = mArg.msize();
1595  tint nOutInfo = 0;
1596 
1597  basic_rvector <float> mD (m);
1598  basic_rvector <float> vOffDiag (_cvm_max<tint>(1, m - 1));
1599  basic_srmatrix <float> mQ (m);
1600  basic_srmatrix <float> mPT (m);
1601  basic_srmatrix <float> mC (1);
1602  basic_rvector <float> vwork (2 * m);
1603  basic_srbmatrix<float> mA (mArg);
1604 
1605  SGBBRD (Chars::pB(),
1606 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1607  1,
1608 #endif
1609  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
1610  mD, vOffDiag,
1611  mQ, mQ._pm(),
1612  mPT, mPT._pm(),
1613  mC, mC._pm(),
1614  vwork, &nOutInfo);
1615  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1616 
1617  basic_rvector<float> vwork2 (m * 4);
1618  SBDSQR (Chars::pU(),
1619 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1620  1,
1621 #endif
1622  &m,
1623  &m,
1624  &m,
1625  &zero,
1626  mD, vOffDiag,
1627  mPT, mPT._pm(),
1628  mQ, mQ._pm(),
1629  mC, mC._pm(),
1630  vwork2, &nOutInfo);
1631  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1632  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "SBDSQR", __FILE__, __LINE__);
1633 
1634  for (tint i = CVM0; i < m + CVM0; ++i) {
1635  if (mD[i] > threshold) {
1636  mPT[i] /= mD[i]; // multiplying V by S^-1
1637  }
1638  else {
1639  mPT[i].set(0.F);
1640  }
1641  }
1642  mX.mult (~mPT, ~mQ);
1643 }
1644 
1645 template<>
1646 CVM_API void
1647 __pinv<double, basic_srbmatrix<double>, basic_rmatrix<double> >
1648  (basic_rmatrix<double>& mX,
1649  const basic_srbmatrix<double>& mArg, double threshold) throw (cvmexception)
1650 {
1651  const tint zero(0);
1652  const tint m = mArg.msize();
1653  tint nOutInfo = 0;
1654 
1655  basic_rvector <double> mD (m);
1656  basic_rvector <double> vOffDiag (_cvm_max<tint>(1, m - 1));
1657  basic_srmatrix <double> mQ (m);
1658  basic_srmatrix <double> mPT (m);
1659  basic_srmatrix <double> mC (1);
1660  basic_rvector <double> vwork (2 * m);
1661  basic_srbmatrix<double> mA (mArg);
1662 
1663  DGBBRD (Chars::pB(),
1664 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1665  1,
1666 #endif
1667  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
1668  mD, vOffDiag,
1669  mQ, mQ._pm(),
1670  mPT, mPT._pm(),
1671  mC, mC._pm(),
1672  vwork, &nOutInfo);
1673  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1674 
1675  basic_rvector<double> vwork2 (m * 4);
1676  DBDSQR (Chars::pU(),
1677 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1678  1,
1679 #endif
1680  &m,
1681  &m,
1682  &m,
1683  &zero,
1684  mD, vOffDiag,
1685  mPT, mPT._pm(),
1686  mQ, mQ._pm(),
1687  mC, mC._pm(),
1688  vwork2, &nOutInfo);
1689  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1690  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "DBDSQR", __FILE__, __LINE__);
1691 
1692  for (tint i = CVM0; i < m + CVM0; ++i) {
1693  if (mD[i] > threshold) {
1694  mPT[i] /= mD[i]; // multiplying V by S^-1
1695  }
1696  else {
1697  mPT[i].set(0.);
1698  }
1699  }
1700  mX.mult (~mPT, ~mQ);
1701 }
1702 
1703 
1705 template<>
1706 CVM_API void
1707 __pinv<float, basic_scbmatrix<float, std::complex<float> >, basic_cmatrix<float, std::complex<float> > >
1708  (basic_cmatrix<float, std::complex<float> >& mX,
1709  const basic_scbmatrix<float, std::complex<float> >& mArg, float threshold) throw (cvmexception)
1710 {
1711  const tint zero(0);
1712  const tint m = mArg.msize();
1713  tint nOutInfo = 0;
1714 
1715  basic_rvector <float> mD (m);
1716  basic_rvector <float> vOffDiag (_cvm_max<tint>(1, m - 1));
1717  basic_scmatrix <float, std::complex<float> > mQ (m);
1718  basic_scmatrix <float, std::complex<float> > mPT (m);
1719  basic_scmatrix <float, std::complex<float> > mC (1);
1720  basic_cvector <float, std::complex<float> > vwork (2 * m);
1721  basic_scbmatrix<float, std::complex<float> > mA (mArg);
1722  basic_rvector <float> vRwork (m);
1723 
1724  CGBBRD (Chars::pB(),
1725 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1726  1,
1727 #endif
1728  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
1729  mD, vOffDiag,
1730  mQ, mQ._pm(),
1731  mPT, mPT._pm(),
1732  mC, mC._pm(),
1733  vwork, vRwork, &nOutInfo);
1734  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1735 
1736  basic_rvector<float> vwork2 (m * 4);
1737  CBDSQR (Chars::pU(),
1738 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1739  1,
1740 #endif
1741  &m,
1742  &m,
1743  &m,
1744  &zero,
1745  mD, vOffDiag,
1746  mPT, mPT._pm(),
1747  mQ, mQ._pm(),
1748  mC, mC._pm(),
1749  vwork2, &nOutInfo);
1750  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1751  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "CBDSQR", __FILE__, __LINE__);
1752 
1753  for (tint i = CVM0; i < m + CVM0; ++i) {
1754  if (mD[i] > threshold) {
1755  mPT[i] /= mD[i]; // multiplying V by S^-1
1756  }
1757  else {
1758  mPT[i].set(0.F);
1759  }
1760  }
1761  mX.mult (~mPT, ~mQ);
1762 }
1763 
1764 template<>
1765 CVM_API void
1766 __pinv<double, basic_scbmatrix<double, std::complex<double> >, basic_cmatrix<double, std::complex<double> > >
1767  (basic_cmatrix<double, std::complex<double> >& mX,
1768  const basic_scbmatrix<double, std::complex<double> >& mArg, double threshold) throw (cvmexception)
1769 {
1770  const tint zero(0);
1771  const tint m = mArg.msize();
1772  tint nOutInfo = 0;
1773 
1774  basic_rvector <double> mD (m);
1775  basic_rvector <double> vOffDiag (_cvm_max<tint>(1, m - 1));
1776  basic_scmatrix <double, std::complex<double> > mQ (m);
1777  basic_scmatrix <double, std::complex<double> > mPT (m);
1778  basic_scmatrix <double, std::complex<double> > mC (1);
1779  basic_cvector <double, std::complex<double> > vwork (2 * m);
1780  basic_rvector <double> vRwork (m);
1781  basic_scbmatrix<double, std::complex<double> > mA (mArg);
1782 
1783  ZGBBRD (Chars::pB(),
1784 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1785  1,
1786 #endif
1787  mA._pm(), mA._pn(), &zero, mA._pl(), mA._pu(), mA, mA._pld(),
1788  mD, vOffDiag,
1789  mQ, mQ._pm(),
1790  mPT, mPT._pm(),
1791  mC, mC._pm(),
1792  vwork, vRwork, &nOutInfo);
1793  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1794 
1795  basic_rvector<double> vwork2 (m * 4);
1796  ZBDSQR (Chars::pU(),
1797 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
1798  1,
1799 #endif
1800  &m,
1801  &m,
1802  &m,
1803  &zero,
1804  mD, vOffDiag,
1805  mPT, mPT._pm(),
1806  mQ, mQ._pm(),
1807  mC, mC._pm(),
1808  vwork2, &nOutInfo);
1809  _check_negative(CVM_WRONGMKLARG, nOutInfo);
1810  _check_positive(CVM_CONVERGENCE_ERROR, nOutInfo, "ZBDSQR", __FILE__, __LINE__);
1811 
1812  for (tint i = CVM0; i < m + CVM0; ++i) {
1813  if (mD[i] > threshold) {
1814  mPT[i] /= mD[i]; // multiplying V by S^-1
1815  }
1816  else {
1817  mPT[i].set(0.);
1818  }
1819  }
1820  mX.mult (~mPT, ~mQ);
1821 }
1823 
1824 CVM_NAMESPACE_END
1825 
Generated on Tue Jul 22 2014 18:52:09 for CVM Class Library by   doxygen 1.8.1.2