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
scmatrix.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 void
16 __exp<basic_scmatrix<float, std::complex<float> >, float>
17  (basic_scmatrix<float, std::complex<float> >& m,
18  const basic_scmatrix<float, std::complex<float> >& mArg,
19  float tol) throw (cvmexception)
20 {
21  tint nR = 0, nI = 0, nQ = 0, nJ = 0;
22  const tint mnM = m.msize();
23  _check_ne(CVM_SIZESMISMATCH, mnM, mArg.msize());
24 
25  basic_scmatrix<float, std::complex<float> > mTmp;
26  const std::complex<float>* pd = mArg._pd();
27 
28  if (pd == m.get())
29  {
30  mTmp << mArg;
31  pd = mTmp;
32  }
33 
34  CMEXPC (&mnM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), &tol, &nR, &nI, &nQ, &nJ);
35  basic_cvector<float, std::complex<float> > vR (nR);
36  basic_array<tint,tint> vI (nI);
37 
38  const tint issymm = 0;
39  std::complex<float> work_dummy(0.F);
40  const tint lwork_dummy = 0;
41 
42  CMEXP (&mnM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), m, m._pld(),
43  vR, vI, &nR, &nI, &nQ, &nJ, &issymm, &work_dummy, &lwork_dummy);
44 }
45 
46 template<>
47 CVM_API void
48 __exp<basic_scmatrix<double, std::complex<double> >, double>
49  (basic_scmatrix<double, std::complex<double> >& m,
50  const basic_scmatrix<double, std::complex<double> >& mArg,
51  double tol) throw (cvmexception)
52 {
53  tint nR = 0, nI = 0, nQ = 0, nJ = 0;
54  const tint mnM = m.msize();
55  _check_ne(CVM_SIZESMISMATCH, mnM, mArg.msize());
56 
57  basic_scmatrix<double, std::complex<double> > mTmp;
58  const std::complex<double>* pd = mArg._pd();
59 
60  if (pd == m.get())
61  {
62  mTmp << mArg;
63  pd = mTmp;
64  }
65 
66  ZMEXPC (&mnM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), &tol, &nR, &nI, &nQ, &nJ);
67  basic_cvector<double, std::complex<double> > vR (nR);
68  basic_array<tint,tint> vI (nI);
69 
70  const tint issymm = 0;
71  std::complex<double> work_dummy(0.);
72  const tint lwork_dummy = 0;
73 
74  ZMEXP (&mnM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), m, m._pld(),
75  vR, vI, &nR, &nI, &nQ, &nJ, &issymm, &work_dummy, &lwork_dummy);
76 }
77 
78 template<>
79 CVM_API void
80 __exp_symm<basic_schmatrix<float, std::complex<float> >, float>
81  (basic_schmatrix<float, std::complex<float> >& m,
82  const basic_schmatrix<float, std::complex<float> >& mArg,
83  float tol) throw (cvmexception)
84 {
85  tint nR = 0, nI = 0, nQ = 0, nJ = 0;
86  const tint nM = m.msize();
87  _check_ne(CVM_SIZESMISMATCH, nM, mArg.msize());
88 
89  basic_schmatrix<float, std::complex<float> > mTmp;
90  const std::complex<float>* pd = mArg._pd();
91 
92  if (pd == m.get())
93  {
94  mTmp << mArg;
95  pd = mTmp;
96  }
97 
98  CMEXPC (&nM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), &tol, &nR, &nI, &nQ, &nJ);
99  basic_cvector<float, std::complex<float> > vR (nR);
100  basic_array<tint,tint> vI (nI);
101 
102  const tint ishem = 1;
103  const tint lwork = 64 * nM;
104  basic_cvector<float, std::complex<float> > work (lwork);
105 
106  CMEXP (&nM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), m, m._pld(),
107  vR, vI, &nR, &nI, &nQ, &nJ, &ishem, work, &lwork);
108 }
109 
110 template<>
111 CVM_API void
112 __exp_symm<basic_schmatrix<double, std::complex<double> >, double>
113  (basic_schmatrix<double, std::complex<double> >& m,
114  const basic_schmatrix<double, std::complex<double> >& mArg,
115  double tol) throw (cvmexception)
116 {
117  tint nR = 0, nI = 0, nQ = 0, nJ = 0;
118  const tint nM = m.msize();
119  _check_ne(CVM_SIZESMISMATCH, nM, mArg.msize());
120 
121  basic_schmatrix<double, std::complex<double> > mTmp;
122  const std::complex<double>* pd = mArg._pd();
123 
124  if (pd == m.get())
125  {
126  mTmp << mArg;
127  pd = mTmp;
128  }
129 
130  ZMEXPC (&nM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), &tol, &nR, &nI, &nQ, &nJ);
131  basic_cvector<double, std::complex<double> > vR (nR);
132  basic_array<tint,tint> vI (nI);
133 
134  const tint ishem = 1;
135  const tint lwork = 64 * nM;
136  basic_cvector<double, std::complex<double> > work (lwork);
137 
138  ZMEXP (&nM, pd, pd == m.get() ? mTmp._pld() : mArg._pldm(), m, m._pld(),
139  vR, vI, &nR, &nI, &nQ, &nJ, &ishem, work, &lwork);
140 }
141 
142 template<>
143 CVM_API void
144 __cond_num<float, basic_scmatrix<float, std::complex<float> > >
145  (const basic_scmatrix<float, std::complex<float> >& mArg, float& dCond) throw (cvmexception)
146 {
147  dCond = 0.F;
148  const tint mnM = mArg.msize();
149  tint nOutInfo = 0;
150  basic_scmatrix<float, std::complex<float> > mA (mArg);
151  basic_cvector<float, std::complex<float> > work (mnM * 2);
152  basic_rvector<float> rwork (mnM * 2);
153  basic_array<tint,tint> iwork (mnM);
154 
155  const float rNorm = mA.norminf();
156  CGETRF (&mnM, &mnM, mA, mA._pld(), iwork, &nOutInfo);
157 
158  _check_negative(CVM_WRONGMKLARG, nOutInfo);
159  if (nOutInfo == 0)
160  {
161 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
162  CGECON (Chars::pI(), 1,
163 #else
164  CGECON (Chars::pI(),
165 #endif
166  &mnM, mA, mA._pld(), &rNorm, &dCond, work, rwork, &nOutInfo);
167  }
168 }
169 
170 template<>
171 CVM_API void
172 __cond_num<double, basic_scmatrix<double, std::complex<double> > >
173  (const basic_scmatrix<double, std::complex<double> >& mArg, double& dCond) throw (cvmexception)
174 {
175  dCond = 0.;
176  const tint mnM = mArg.msize();
177  tint nOutInfo = 0;
178  basic_scmatrix<double, std::complex<double> > mA (mArg);
179  basic_cvector<double, std::complex<double> > work (mnM * 2);
180  basic_rvector<double> rwork (mnM * 2);
181  basic_array<tint,tint> iwork (mnM);
182 
183  const double rNorm = mA.norminf();
184  ZGETRF (&mnM, &mnM, mA, mA._pld(), iwork, &nOutInfo);
185 
186  _check_negative(CVM_WRONGMKLARG, nOutInfo);
187  if (nOutInfo == 0)
188  {
189 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
190  ZGECON (Chars::pI(), 1,
191 #else
192  ZGECON (Chars::pI(),
193 #endif
194  &mnM, mA, mA._pld(), &rNorm, &dCond, work, rwork, &nOutInfo);
195  }
196 }
197 
198 template <>
199 CVM_API void
200 __inv<basic_scmatrix<float, std::complex<float> > >
201  (basic_scmatrix<float, std::complex<float> >& m,
202  const basic_scmatrix<float, std::complex<float> >& mArg) throw (cvmexception)
203 {
204  const tint mnM = m.msize();
205  _check_ne(CVM_SIZESMISMATCH, mnM, mArg.msize());
206 
207  if (mnM == 1)
208  {
209  if (_abs (mArg(CVM0,CVM0)) <= basic_cvmMachMin<float>()) {
210  throw cvmexception (CVM_SINGULARMATRIX, 1);
211  }
212  m(CVM0,CVM0) = 1.F / mArg(CVM0,CVM0);
213  }
214  else
215  {
216  basic_array<tint,tint> nPivots (mnM);
217  m.low_up (mArg, nPivots);
218 
219  tint lWork = -1;
220  tint nOutInfo = 0;
221  std::complex<float> dWork;
222  CGETRI (&mnM, m, m._pld(), nPivots, &dWork, &lWork, &nOutInfo);
223  _check_negative(CVM_WRONGMKLARG, nOutInfo);
224  lWork = static_cast<tint>(dWork.real());
225  basic_cvector<float, std::complex<float> > vWork (lWork);
226 
227  CGETRI (&mnM, m, m._pld(), nPivots, vWork, &lWork, &nOutInfo);
228  _check_negative(CVM_WRONGMKLARG, nOutInfo);
229  _check_positive(CVM_SINGULARMATRIX, nOutInfo);
230  }
231 }
232 
233 template<>
234 CVM_API void
235 __inv<basic_scmatrix<double, std::complex<double> > >
236  (basic_scmatrix<double, std::complex<double> >& m,
237  const basic_scmatrix<double, std::complex<double> >& mArg) throw (cvmexception)
238 {
239  const tint mnM = m.msize();
240  _check_ne(CVM_SIZESMISMATCH, mnM, mArg.msize());
241 
242  if (mnM == 1)
243  {
244  if (_abs (mArg(CVM0,CVM0)) <= basic_cvmMachMin<double>()) {
245  throw cvmexception (CVM_SINGULARMATRIX, 1);
246  }
247  m(CVM0,CVM0) = 1. / mArg(CVM0,CVM0);
248  }
249  else
250  {
251  basic_array<tint,tint> nPivots (mnM);
252  m.low_up (mArg, nPivots);
253 
254  tint lWork = -1;
255  tint nOutInfo = 0;
256  std::complex<double> dWork;
257  ZGETRI (&mnM, m, m._pld(), nPivots, &dWork, &lWork, &nOutInfo);
258  lWork = static_cast<tint>(dWork.real());
259  basic_cvector<double, std::complex<double> > vWork (lWork);
260 
261  ZGETRI (&mnM, m, m._pld(), nPivots, vWork, &lWork, &nOutInfo);
262  _check_negative(CVM_WRONGMKLARG, nOutInfo);
263  _check_positive(CVM_SINGULARMATRIX, nOutInfo);
264  }
265 }
266 
267 template<>
268 CVM_API void
269 __inv<basic_schmatrix<float, std::complex<float> > >
270  (basic_schmatrix<float, std::complex<float> >& m,
271  const basic_schmatrix<float, std::complex<float> >& mArg) throw (cvmexception)
272 {
273  const tint nM = m.msize();
274  _check_ne(CVM_SIZESMISMATCH, nM, mArg.msize());
275 
276  if (nM == 1)
277  {
278  if (_abs (mArg(CVM0,CVM0)) <= basic_cvmMachMin<float>()) {
279  throw cvmexception (CVM_SINGULARMATRIX, 1);
280  }
281  m(CVM0,CVM0) = 1.F / mArg(CVM0,CVM0);
282  }
283  else
284  {
285  bool bPositiveDefinite = false;
286  tint nOutInfo = 0;
287  basic_array<tint,tint> nPivots (nM);
288 
289  m._factorize (mArg, nPivots, bPositiveDefinite);
290 
291  if (bPositiveDefinite)
292  {
293  CPOTRI (Chars::pU(),
294 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
295  1,
296 #endif
297  &nM, m, m._pld(), &nOutInfo);
298  _check_negative(CVM_WRONGMKLARG, nOutInfo);
299  _check_positive(CVM_WRONGCHOLESKYFACTOR, nOutInfo);
300  }
301  else
302  {
303  basic_cvector<float, std::complex<float> > vWork (nM);
304  CHETRI (Chars::pU(),
305 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
306  1,
307 #endif
308  &nM, m, m._pld(), nPivots, vWork, &nOutInfo);
309  _check_negative(CVM_WRONGMKLARG, nOutInfo);
310  _check_positive(CVM_WRONGBUNCHKAUFMANFACTOR, nOutInfo);
311  }
312  m._flip();
313  }
314 }
315 
316 template<>
317 CVM_API void
318 __inv<basic_schmatrix<double, std::complex<double> > >
319  (basic_schmatrix<double, std::complex<double> >& m,
320  const basic_schmatrix<double, std::complex<double> >& mArg) throw (cvmexception)
321 {
322  const tint nM = m.msize();
323  _check_ne(CVM_SIZESMISMATCH, nM, mArg.msize());
324 
325  if (nM == 1)
326  {
327  if (_abs (mArg(CVM0,CVM0)) <= basic_cvmMachMin<double>()) {
328  throw cvmexception (CVM_SINGULARMATRIX, 1);
329  }
330  m(CVM0,CVM0) = 1. / mArg(CVM0,CVM0);
331  }
332  else
333  {
334  bool bPositiveDefinite = false;
335  tint nOutInfo = 0;
336  basic_array<tint,tint> nPivots (nM);
337 
338  m._factorize (mArg, nPivots, bPositiveDefinite);
339 
340  if (bPositiveDefinite)
341  {
342  ZPOTRI (Chars::pU(),
343 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
344  1,
345 #endif
346  &nM, m, m._pld(), &nOutInfo);
347  _check_negative(CVM_WRONGMKLARG, nOutInfo);
348  _check_positive(CVM_WRONGCHOLESKYFACTOR, nOutInfo);
349  }
350  else
351  {
352  basic_cvector<double, std::complex<double> > vWork (nM);
353  ZHETRI (Chars::pU(),
354 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
355  1,
356 #endif
357  &nM, m, m._pld(), nPivots, vWork, &nOutInfo);
358  _check_negative(CVM_WRONGMKLARG, nOutInfo);
359  _check_positive(CVM_WRONGBUNCHKAUFMANFACTOR, nOutInfo);
360  }
361  m._flip();
362  }
363 }
364 
366 template<>
367 CVM_API void
368 __polynom<std::complex<float>, basic_cvector<float, std::complex<float> > >
369  (std::complex<float>* mpd, tint ldP,
370  tint mnM,
371  const std::complex<float>* pd, tint ldA,
372  const basic_cvector<float, std::complex<float> >& v)
373 {
374  basic_cvector<float, std::complex<float> > vWork (NPOLY (&mnM, v._psize()));
375  CPOLY (&mnM, pd, &ldA, v._psize(), v, mpd, &ldP, vWork);
376 }
377 
378 template<>
379 CVM_API void
380 __polynom<std::complex<double>, basic_cvector<double, std::complex<double> > >
381  (std::complex<double>* mpd, tint ldP,
382  tint mnM,
383  const std::complex<double>* pd, tint ldA,
384  const basic_cvector<double, std::complex<double> >& v)
385 {
386  basic_cvector<double, std::complex<double> > vWork (NPOLY (&mnM, v._psize()));
387  ZPOLY (&mnM, pd, &ldA, v._psize(), v, mpd, &ldP, vWork);
388 }
389 
390 // internal solver
391 // do not forget to make pX equal to pB before call
392 // in case of vectors passed increment MUST be 1 (since those methods assume matrices)
393 template<>
394 CVM_API void
395 __solve<float, std::complex<float>, basic_scmatrix<float, std::complex<float> > >
396  (const basic_scmatrix<float, std::complex<float> >& m,
397  tint nrhs,
398  const std::complex<float>* pB, tint ldB,
399  std::complex<float>* pX, tint ldX,
400  float& dErr,
401  const std::complex<float>* pLU, const tint* pPivots, int transp_mode) throw (cvmexception)
402 {
403  const tint mnM = m.msize();
404  const bool bGivenLU = pLU != nullptr && pPivots != nullptr;
405  tint nOutInfo = 0;
406  basic_rvector<float> vFerr (nrhs);
407  basic_rvector<float> vBerr (nrhs);
408  basic_cvector<float, std::complex<float> > vWork (2 * mnM);
409  basic_rvector<float> rWork (mnM);
410  basic_array<tint,tint> nPivots (mnM);
411  // 0 - none, 1 - just transposed, 2 - conjugated
412  const char* transp = transp_mode == 0 ? Chars::pN() : (transp_mode == 1 ? Chars::pT() : Chars::pC());
413 
414  if (bGivenLU) nPivots.assign (pPivots);
415  basic_scmatrix<float, std::complex<float> > mLU (mnM);
416  if (bGivenLU)
417  {
418  mLU.assign (pLU);
419  }
420  else
421  {
422  mLU = m.low_up (nPivots);
423  }
424  dErr = 0.F;
425 
426  CGETRS (transp,
427 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
428  1,
429 #endif
430  &mnM, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
431 
432  _check_negative(CVM_WRONGMKLARG, nOutInfo);
433 
434  CGERFS (transp,
435 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
436  1,
437 #endif
438  &mnM, &nrhs,
439  m, m._pld(),
440  mLU, mLU._pld(),
441  nPivots,
442  pB, &ldB,
443  pX, &ldX,
444  vFerr, vBerr, vWork, rWork, &nOutInfo);
445 
446  _check_negative(CVM_WRONGMKLARG, nOutInfo);
447  dErr = vFerr.norminf();
448 }
449 
450 template<>
451 CVM_API void
452 __solve<double, std::complex<double>, basic_scmatrix<double, std::complex<double> > >
453  (const basic_scmatrix<double, std::complex<double> >& m,
454  tint nrhs,
455  const std::complex<double>* pB, tint ldB,
456  std::complex<double>* pX, tint ldX,
457  double& dErr,
458  const std::complex<double>* pLU, const tint* pPivots, int transp_mode) throw (cvmexception)
459 {
460  const tint mnM = m.msize();
461  const bool bGivenLU = pLU != nullptr && pPivots != nullptr;
462  tint nOutInfo = 0;
463  basic_rvector<double> vBerr (nrhs);
464  basic_rvector<double> vFerr (nrhs);
465  basic_cvector<double, std::complex<double> > vWork (2 * mnM);
466  basic_rvector<double> rWork (mnM);
467  basic_array<tint,tint> nPivots (mnM);
468  // 0 - none, 1 - just transposed, 2 - conjugated
469  const char* transp = transp_mode == 0 ? Chars::pN() : (transp_mode == 1 ? Chars::pT() : Chars::pC());
470 
471  if (bGivenLU) nPivots.assign (pPivots);
472  basic_scmatrix<double, std::complex<double> > mLU (mnM);
473  if (bGivenLU)
474  {
475  mLU.assign (pLU);
476  }
477  else
478  {
479  mLU = m.low_up (nPivots);
480  }
481  dErr = 0.;
482 
483  ZGETRS (transp,
484 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
485  1,
486 #endif
487  &mnM, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
488 
489  _check_negative(CVM_WRONGMKLARG, nOutInfo);
490 
491  ZGERFS (transp,
492 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
493  1,
494 #endif
495  &mnM, &nrhs,
496  m, m._pld(),
497  mLU, mLU._pld(),
498  nPivots,
499  pB, &ldB,
500  pX, &ldX,
501  vFerr, vBerr, vWork, rWork, &nOutInfo);
502 
503  _check_negative(CVM_WRONGMKLARG, nOutInfo);
504  dErr = vFerr.norminf();
505 }
506 
507 template<>
508 CVM_API void
509 __solve<float, std::complex<float>, basic_scbmatrix<float, std::complex<float> > >
510  (const basic_scbmatrix<float, std::complex<float> >& m,
511  tint nrhs,
512  const std::complex<float>* pB, tint ldB,
513  std::complex<float>* pX, tint ldX,
514  float& dErr,
515  const std::complex<float>* pLU, const tint* pPivots, int transp_mode) throw (cvmexception)
516 {
517  const tint mnM = m.msize();
518  const tint mnKL = m.lsize();
519  const tint mnKU = m.usize();
520  const bool bGivenLU = pLU != nullptr && pPivots != nullptr;
521  tint nOutInfo = 0;
522  basic_rvector<float> vFerr (nrhs);
523  basic_rvector<float> vBerr (nrhs);
524  basic_cvector<float, std::complex<float> > vWork (2 * mnM);
525  basic_rvector<float> rWork (mnM);
526  basic_array<tint,tint>nPivots (mnM);
527  // 0 - none, 1 - just transposed, 2 - conjugated
528  const char* transp = transp_mode == 0 ? Chars::pN() : (transp_mode == 1 ? Chars::pT() : Chars::pC());
529 
530  if (bGivenLU) nPivots.assign (pPivots);
531  basic_scbmatrix<float, std::complex<float> > mLU (mnM, mnKL, mnKL + mnKU);
532  if (bGivenLU)
533  {
534  mLU.assign (pLU);
535  }
536  else
537  {
538  mLU = m.low_up (nPivots);
539  }
540  dErr = 0.F;
541 
542  CGBTRS (transp,
543 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
544  1,
545 #endif
546  &mnM, &mnKL, &mnKU, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
547 
548  _check_negative(CVM_WRONGMKLARG, nOutInfo);
549 
550  CGBRFS (transp,
551 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
552  1,
553 #endif
554  &mnM, &mnKL, &mnKU, &nrhs,
555  m, m._pld(),
556  mLU, mLU._pld(),
557  nPivots,
558  pB, &ldB,
559  pX, &ldX,
560  vFerr, vBerr, vWork, rWork, &nOutInfo);
561 
562  _check_negative(CVM_WRONGMKLARG, nOutInfo);
563  dErr = vFerr.norminf();
564 }
565 
566 template<>
567 CVM_API void
568 __solve<double, std::complex<double>, basic_scbmatrix<double, std::complex<double> > >
569  (const basic_scbmatrix<double, std::complex<double> >& m,
570  tint nrhs,
571  const std::complex<double>* pB, tint ldB,
572  std::complex<double>* pX, tint ldX,
573  double& dErr,
574  const std::complex<double>* pLU, const tint* pPivots, int transp_mode) throw (cvmexception)
575 {
576  const tint mnM = m.msize();
577  const tint mnKL = m.lsize();
578  const tint mnKU = m.usize();
579  const bool bGivenLU = pLU != nullptr && pPivots != nullptr;
580  tint nOutInfo = 0;
581  basic_rvector<double> vFerr (nrhs);
582  basic_rvector<double> vBerr (nrhs);
583  basic_cvector<double, std::complex<double> > vWork (2 * mnM);
584  basic_rvector<double> rWork (mnM);
585  basic_array<tint,tint> nPivots (mnM);
586  // 0 - none, 1 - just transposed, 2 - conjugated
587  const char* transp = transp_mode == 0 ? Chars::pN() : (transp_mode == 1 ? Chars::pT() : Chars::pC());
588 
589  if (bGivenLU) nPivots.assign (pPivots);
590  basic_scbmatrix<double, std::complex<double> > mLU (mnM, mnKL, mnKL + mnKU);
591  if (bGivenLU)
592  {
593  mLU.assign (pLU);
594  }
595  else
596  {
597  mLU = m.low_up (nPivots);
598  }
599  dErr = 0.;
600 
601  ZGBTRS (transp,
602 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
603  1,
604 #endif
605  &mnM, &mnKL, &mnKU, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
606 
607  _check_negative(CVM_WRONGMKLARG, nOutInfo);
608 
609  ZGBRFS (transp,
610 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
611  1,
612 #endif
613  &mnM, &mnKL, &mnKU, &nrhs,
614  m, m._pld(),
615  mLU, mLU._pld(),
616  nPivots,
617  pB, &ldB,
618  pX, &ldB,
619  vFerr, vBerr, vWork, rWork, &nOutInfo);
620 
621  _check_negative(CVM_WRONGMKLARG, nOutInfo);
622  dErr = vFerr.norminf();
623 }
624 
625 
626 // internal solver
627 // don't forget to make pX equal to pB before call
628 template<>
629 CVM_API void
630 __solve<float, std::complex<float>, basic_schmatrix<float, std::complex<float> > >
631  (const basic_schmatrix<float, std::complex<float> >& m,
632  tint nrhs,
633  const std::complex<float>* pB, tint ldB,
634  std::complex<float>* pX, tint ldX,
635  float& dErr,
636  const std::complex<float>* pLU, const tint* pPivots, int) throw (cvmexception)
637 {
638  const tint nM = m.msize();
639  const bool bCholeskyGiven = pLU != nullptr && pPivots == nullptr; // no pivots means Cholesky
640  const bool bBunchKaufmanGiven = pLU != nullptr && pPivots != nullptr;
641  const bool bCalc = !bCholeskyGiven && !bBunchKaufmanGiven;
642  bool bPositiveDefinite = bCholeskyGiven;
643 
644  tint nOutInfo = 0;
645  basic_rvector<float> vBerr (nrhs);
646  basic_rvector<float> vFerr (nrhs);
647  basic_cvector<float, std::complex<float> > vWork (2 * nM);
648  basic_rvector<float> vrWork (nM);
649  basic_array<tint,tint> nPivots (nM);
650 
651  if (bBunchKaufmanGiven) nPivots.assign (pPivots);
652  basic_schmatrix<float, std::complex<float> > mLU(nM);
653  if (bCalc)
654  {
655  mLU._factorize (m, nPivots, bPositiveDefinite);
656  }
657  else
658  {
659  mLU.assign (pLU);
660  }
661  dErr = 0.F;
662 
663  if (bPositiveDefinite)
664  {
665  CPOTRS (Chars::pU(),
666 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
667  1,
668 #endif
669  &nM, &nrhs, mLU, mLU._pld(), pX, &ldX, &nOutInfo);
670  _check_negative(CVM_WRONGMKLARG, nOutInfo);
671 
672  CPORFS (Chars::pU(),
673 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
674  1,
675 #endif
676  &nM, &nrhs,
677  m, m._pld(),
678  mLU, mLU._pld(),
679  pB, &ldB,
680  pX, &ldX,
681  vFerr, vBerr,
682  vWork, vrWork, &nOutInfo);
683  _check_negative(CVM_WRONGMKLARG, nOutInfo);
684  }
685  else
686  {
687  CHETRS (Chars::pU(),
688 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
689  1,
690 #endif
691  &nM, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
692  _check_negative(CVM_WRONGMKLARG, nOutInfo);
693 
694  CHERFS (Chars::pU(),
695 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
696  1,
697 #endif
698  &nM, &nrhs,
699  m, m._pld(),
700  mLU, mLU._pld(),
701  nPivots,
702  pB, &ldB,
703  pX, &ldX,
704  vFerr, vBerr,
705  vWork, vrWork, &nOutInfo);
706  _check_negative(CVM_WRONGMKLARG, nOutInfo);
707  }
708 
709  dErr = vFerr.norminf();
710 }
711 
712 template<>
713 CVM_API void
714 __solve<double, std::complex<double>, basic_schmatrix<double, std::complex<double> > >
715  (const basic_schmatrix<double, std::complex<double> >& m,
716  tint nrhs,
717  const std::complex<double>* pB, tint ldB,
718  std::complex<double>* pX, tint ldX,
719  double& dErr,
720  const std::complex<double>* pLU, const tint* pPivots, int) throw (cvmexception)
721 {
722  const tint nM = m.msize();
723  const bool bCholeskyGiven = pLU != nullptr && pPivots == nullptr; // no pivots means Cholesky
724  const bool bBunchKaufmanGiven = pLU != nullptr && pPivots != nullptr;
725  const bool bCalc = !bCholeskyGiven && !bBunchKaufmanGiven;
726  bool bPositiveDefinite = bCholeskyGiven;
727 
728  tint nOutInfo = 0;
729  basic_rvector<double> vBerr (nrhs);
730  basic_rvector<double> vFerr (nrhs);
731  basic_cvector<double, std::complex<double> > vWork (2 * nM);
732  basic_rvector<double> vrWork (nM);
733  basic_array<tint,tint> nPivots (nM);
734 
735  if (bBunchKaufmanGiven) nPivots.assign (pPivots);
736  basic_schmatrix<double, std::complex<double> > mLU(nM);
737  if (bCalc)
738  {
739  mLU._factorize (m, nPivots, bPositiveDefinite);
740  }
741  else
742  {
743  mLU.assign (pLU);
744  }
745  dErr = 0.;
746 
747  if (bPositiveDefinite)
748  {
749  ZPOTRS (Chars::pU(),
750 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
751  1,
752 #endif
753  &nM, &nrhs, mLU, mLU._pld(), pX, &ldX, &nOutInfo);
754  _check_negative(CVM_WRONGMKLARG, nOutInfo);
755 
756  ZPORFS (Chars::pU(),
757 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
758  1,
759 #endif
760  &nM, &nrhs,
761  m, m._pld(),
762  mLU, mLU._pld(),
763  pB, &ldB,
764  pX, &ldX,
765  vFerr, vBerr,
766  vWork, vrWork, &nOutInfo);
767  _check_negative(CVM_WRONGMKLARG, nOutInfo);
768  }
769  else
770  {
771  ZHETRS (Chars::pU(),
772 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
773  1,
774 #endif
775  &nM, &nrhs, mLU, mLU._pld(), nPivots, pX, &ldX, &nOutInfo);
776  _check_negative(CVM_WRONGMKLARG, nOutInfo);
777 
778  ZHERFS (Chars::pU(),
779 #if defined (CVM_PASS_STRING_LENGTH_TO_FTN_SUBROUTINES)
780  1,
781 #endif
782  &nM, &nrhs,
783  m, m._pld(),
784  mLU, mLU._pld(),
785  nPivots,
786  pB, &ldB,
787  pX, &ldX,
788  vFerr, vBerr,
789  vWork, vrWork, &nOutInfo);
790  _check_negative(CVM_WRONGMKLARG, nOutInfo);
791  }
792 
793  dErr = vFerr.norminf();
794 }
796 
797 CVM_NAMESPACE_END
798 
Generated on Tue Jul 22 2014 18:52:09 for CVM Class Library by   doxygen 1.8.1.2