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