SLIP  1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
linear_algebra.hpp
Go to the documentation of this file.
1 /*
2  * Copyright(c):
3  * Signal Image and Communications (SIC) Department
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * - University of Poitiers, France http://www.univ-poitiers.fr
6  * - XLIM Institute UMR CNRS 7252 http://www.xlim.fr/
7  *
8  * and
9  *
10  * D2 Fluid, Thermic and Combustion
11  * - University of Poitiers, France http://www.univ-poitiers.fr
12  * - PPRIME Institute - UPR CNRS 3346 http://www.pprime.fr
13  * - ISAE-ENSMA http://www.ensma.fr
14  *
15  * Contributor(s):
16  * The SLIP team,
17  * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
18  * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>,
19  * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>,
20  * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>,
21  * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>,
22  * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
23  *
24  * Description:
25  * The Simple Library of Image Processing (SLIP) is a new image processing
26  * library. It is written in the C++ language following as much as possible
27  * the ISO/ANSI C++ standard. It is consequently compatible with any system
28  * satisfying the ANSI C++ complience. It works on different Unix , Linux ,
29  * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that
30  * was created by the Signal, Image and Communications (SIC) departement of
31  * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic
32  * and Combustion departement of the P', UPR 3346 CNRS Institute of the
33  * University of Poitiers.
34  *
35  * The SLIP Library source code has been registered to the APP (French Agency
36  * for the Protection of Programs) by the University of Poitiers and CNRS,
37  * under registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.
38 
39  * http://www.sic.sp2mi.univ-poitiers.fr/slip/
40  *
41  * This software is governed by the CeCILL-C license under French law and
42  * abiding by the rules of distribution of free software. You can use,
43  * modify and/ or redistribute the software under the terms of the CeCILL-C
44  * license as circulated by CEA, CNRS and INRIA at the following URL
45  * http://www.cecill.info.
46  * As a counterpart to the access to the source code and rights to copy,
47  * modify and redistribute granted by the license, users are provided only
48  * with a limited warranty and the software's author, the holder of the
49  * economic rights, and the successive licensors have only limited
50  * liability.
51  *
52  * In this respect, the user's attention is drawn to the risks associated
53  * with loading, using, modifying and/or developing or reproducing the
54  * software by the user in light of its specific status of free software,
55  * that may mean that it is complicated to manipulate, and that also
56  * therefore means that it is reserved for developers and experienced
57  * professionals having in-depth computer knowledge. Users are therefore
58  * encouraged to load and test the software's suitability as regards their
59  * requirements in conditions enabling the security of their systems and/or
60  * data to be ensured and, more generally, to use and operate it in the
61  * same conditions as regards security.
62  *
63  * The fact that you are presently reading this means that you have had
64  * knowledge of the CeCILL-C license and that you accept its terms.
65  */
66 
67 
74 #ifndef SLIP_LINEAR_ALGEBRA_HPP
75 #define SLIP_LINEAR_ALGEBRA_HPP
76 
77 #include <iterator>
78 #include <algorithm>
79 #include <numeric>
80 #include <cassert>
81 #include <cmath>
82 #include <limits>
83 #include <iostream>
84 #include <complex>
85 #include <stdexcept>
86 #include "complex_cast.hpp"
87 #include "Array2d.hpp"
88 //#include "Vector.hpp"
89 #include "Array.hpp"
90 #include "macros.hpp"
91 #include "arithmetic_op.hpp"
92 #include "linear_algebra_traits.hpp"
93 #include "error.hpp"
94 #include "apply.hpp"
95 #include "norms.hpp"
96 
97 namespace slip
98 {
99 
100 template <class RandomAccessIterator2d,
101  class RandomAccessIterator1,
102  class RandomAccessIterator2>
103 void matrix_vector_multiplies(RandomAccessIterator2d M_upper_left,
104  RandomAccessIterator2d M_bottom_right,
105  RandomAccessIterator1 first1,
106  RandomAccessIterator1 last1,
107  RandomAccessIterator2 result_first,
108  RandomAccessIterator2 result_last);
134 enum DENSE_MATRIX_TYPE
135 {
136  DENSE_MATRIX_FULL,
137  DENSE_MATRIX_DIAGONAL,
138  DENSE_MATRIX_UPPER_BIDIAGONAL,
139  DENSE_MATRIX_LOWER_BIDIAGONAL,
140  DENSE_MATRIX_TRIDIAGONAL,
141  DENSE_MATRIX_UPPER_TRIANGULAR,
142  DENSE_MATRIX_LOWER_TRIANGULAR,
143  DENSE_MATRIX_BANDED,
144  DENSE_MATRIX_UPPER_HESSENBERG,
145  DENSE_MATRIX_LOWER_HESSENBERG,
146  DENSE_MATRIX_NULL
147 };
148 
149 
161  template <typename _F, typename _S, typename _RT, typename _AT>
162  struct saxpy : public std::binary_function<_F, _S, _RT>
163  {
164  saxpy():a_(_AT(1))
165  {}
166  saxpy(const _AT& a):a_(a)
167  {}
168 
169  inline
170  _RT
171  operator()(const _F& __x, const _S& __y) const
172  {
173  return a_ * __x + __y;
174  }
175 
176  _AT a_;
177  };
178 
179 
193  template <class _Tp>
194  struct z1conjz2 : public std::binary_function<_Tp, _Tp, _Tp>
195  {
196  inline
197  _Tp
198  operator()(const _Tp& z1, const _Tp& z2) const
199  {
200  return z1 * slip::conj(z2);
201  }
202  };
203 
217  template <class _Tp>
218  struct conjz1z2 : public std::binary_function<_Tp, _Tp, _Tp>
219  {
220  inline
221  _Tp
222  operator()(const _Tp& z1, const _Tp& z2) const
223  {
224  return slip::conj(z1) * z2;
225  }
226  };
227 
228 
229 
230 
231 
233  /* @{ */
267  template <typename RandomAccessIterator1,
268  typename RandomAccessIterator2>
269  inline
270  typename std::iterator_traits<RandomAccessIterator1>::value_type
271  hermitian_inner_product(RandomAccessIterator1 first1,
272  RandomAccessIterator1 last1,
273  RandomAccessIterator2 first2,
274  typename std::iterator_traits<RandomAccessIterator1>::value_type init)
275  {
276  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
277  return std::inner_product(first1,
278  last1,
279  first2,
280  init,
281  std::plus<value_type>(),
282  slip::conjz1z2<value_type>());
283  }
284 
285 
319  template <typename RandomAccessIterator1,
320  typename RandomAccessIterator2>
321  inline
322  typename std::iterator_traits<RandomAccessIterator1>::value_type
323  inner_product(RandomAccessIterator1 first1,
324  RandomAccessIterator1 last1,
325  RandomAccessIterator2 first2,
326  typename std::iterator_traits<RandomAccessIterator1>::value_type init = typename std::iterator_traits<RandomAccessIterator1>::value_type())
327  {
328  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
329 
330 
331  if(slip::is_complex(value_type()))
332  {
333  return slip::hermitian_inner_product(first1,
334  last1,
335  first2,
336  init);
337  }
338  else
339  {
340  return std::inner_product(first1,
341  last1,
342  first2,
343  init);
344  }
345  }
346 
347 
377  template <typename Vector1,
378  typename Vector2>
379  inline
380  typename Vector1::value_type
381  inner_product(const Vector1& V1,
382  const Vector2& V2)
383  {
384  assert(V1.size() >= V2.size());
385  return slip::inner_product(V1.begin(),V1.end(),V2.begin(),typename Vector1::value_type());
386  }
387 
388 
389 
390 
391 
392 
393 
437  template <typename RandomAccessIterator1,
438  typename RandomAccessIterator2d,
439  typename RandomAccessIterator2>
440  inline
441  typename std::iterator_traits<RandomAccessIterator2>::value_type
442  gen_inner_product(RandomAccessIterator1 first1,
443  RandomAccessIterator1 last1,
444  RandomAccessIterator2d A_upper_left,
445  RandomAccessIterator2d A_bottom_right,
446  RandomAccessIterator2 first2,
447  RandomAccessIterator2 last2,
448  typename std::iterator_traits<RandomAccessIterator2>::value_type init = typename std::iterator_traits<RandomAccessIterator2>::value_type())
449  {
450  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
451  slip::Array<value_type> tmp(last1-first1);
452  slip::matrix_vector_multiplies(A_upper_left,A_bottom_right,first2,last2,tmp.begin(),tmp.end());
453 
454  if(slip::is_complex(value_type()))
455  {
456  return slip::hermitian_inner_product(first1,last1,
457  tmp.begin(),
458  init);
459  }
460  else
461  {
462  return std::inner_product(first1,last1,
463  tmp.begin(),
464  init);
465  }
466  }
467 
468 
504  template <typename Vector1,
505  typename Matrix,
506  typename Vector2>
507  inline
508  typename Vector2::value_type
509  gen_inner_product(const Vector1& X,
510  const Matrix& M,
511  const Vector2& Y)
512 
513  {
514  return slip::gen_inner_product(X.begin(),X.end(),
515  M.upper_left(),M.bottom_right(),
516  Y.begin(),Y.end(),
517  typename Vector2::value_type());
518 
519  }
520 
521 
553  template <typename RandomAccessIterator1,
554  typename RandomAccessIterator2d>
555  inline
556  typename std::iterator_traits<RandomAccessIterator1>::value_type
557  rayleigh(RandomAccessIterator1 first1,
558  RandomAccessIterator1 last1,
559  RandomAccessIterator2d A_upper_left,
560  RandomAccessIterator2d A_bottom_right)
561  {
562 
563  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
564  return slip::gen_inner_product(first1,last1,
565  A_upper_left,A_bottom_right,
566  first1,last1) /
567  slip::inner_product(first1,last1,
568  first1,
569  value_type());
570  }
571 
599  template <typename Vector,
600  typename Matrix>
601  inline
602  typename Vector::value_type
603  rayleigh(const Vector& X,
604  const Matrix& A)
605  {
606  return slip::rayleigh(X.begin(),X.end(),
607  A.upper_left(),A.bottom_right());
608  }
609 
610 /* @} */
612  /* @{ */
613 
628  template <typename AT,
629  typename RandomAccessIterator1,
630  typename RandomAccessIterator2>
631  inline
632  void aXpY(const AT& a,
633  RandomAccessIterator1 first1,
634  RandomAccessIterator1 last1,
635  RandomAccessIterator2 first2)
636  {
637  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type1;
638  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type2;
639  std::transform(first1,last1,
640  first2,
641  first2,
642  slip::saxpy<value_type1,value_type2,value_type2,AT>(a));
643  }
644 
684  template <typename RandomAccessIterator1,
685  typename T,
686  typename RandomAccessIterator2>
687  inline
688  void vector_scalar_multiplies(RandomAccessIterator1 V_first,
689  RandomAccessIterator1 V_last,
690  const T& scal,
691  RandomAccessIterator2 result_first)
692  {
693 
694 
695  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type1;
696 
697  for ( ; V_first != V_last; ++V_first, ++result_first)
698  {
699  *result_first = *V_first * scal;
700  }
701  }
702 
703 
729  template <typename Vector1,
730  typename T,
731  typename Vector2>
732  inline
733  void vector_scalar_multiplies(const Vector1& V,
734  const T& scal,
735  Vector2& Result)
736  {
737  assert(Result.size() <= V.size());
738  slip::vector_scalar_multiplies(V.begin(),V.end(),scal,Result.begin());
739  }
740 
770  template <typename RandomAccessIterator1,
771  typename T,
772  typename RandomAccessIterator2>
773  inline
774  void conj_vector_scalar_multiplies(RandomAccessIterator1 V_first,
775  RandomAccessIterator1 V_last,
776  const T& scal,
777  RandomAccessIterator2 result_first)
778  {
779  for ( ; V_first != V_last; ++V_first, ++result_first)
780  {
781  *result_first = slip::conj(*V_first) * scal;
782  }
783  }
784 
811  template <typename Vector1,
812  typename T,
813  typename Vector2>
814  inline
815  void conj_vector_scalar_multiplies(const Vector1& V,
816  const T& scal,
817  Vector2& Result)
818  {
819  slip::conj_vector_scalar_multiplies(V.begin(),V.end(),scal,Result.begin());
820  }
821 
852  template <class RandomAccessIterator2d1,
853  class T,
854  class RandomAccessIterator2d2>
855  inline
856  void matrix_scalar_multiplies(RandomAccessIterator2d1 M_upper_left,
857  RandomAccessIterator2d1 M_bottom_right,
858  const T& scal,
859  RandomAccessIterator2d2 R_upper_left,
860  RandomAccessIterator2d2 R_bottom_right)
861 
862  {
863 
864  typedef typename RandomAccessIterator2d1::size_type size_type;
865  size_type M_rows = (M_bottom_right - M_upper_left)[0];
866 
867  assert((M_bottom_right - M_upper_left)[0] == (R_bottom_right - R_upper_left)[0]);
868  assert((M_bottom_right - M_upper_left)[1] == (R_bottom_right - R_upper_left)[1]);
869 
870 
871  for(std::size_t i = 0; i < M_rows; ++i)
872  {
873  slip::vector_scalar_multiplies(M_upper_left.row_begin(i),
874  M_upper_left.row_end(i),
875  scal,
876  R_upper_left.row_begin(i));
877  }
878 
879  }
880 
881 
906  template <class Matrix1,
907  class T,
908  class Matrix2>
909  inline
910  void matrix_scalar_multiplies(const Matrix1& M,
911  const T& scal,
912  Matrix2& R)
913 
914  {
915 
916  slip::matrix_scalar_multiplies(M.upper_left(),M.bottom_right(),
917  scal,
918  R.upper_left(),R.bottom_right());
919 
920  }
921 
922 
923 
959  template <class RandomAccessIterator2d,
960  class RandomAccessIterator1,
961  class RandomAccessIterator2>
962  inline
963  void matrix_vector_multiplies(RandomAccessIterator2d M_upper_left,
964  RandomAccessIterator2d M_bottom_right,
965  RandomAccessIterator1 first1,
966  RandomAccessIterator1 last1,
967  RandomAccessIterator2 result_first,
968  RandomAccessIterator2 result_last)
969  {
970  assert((M_bottom_right - M_upper_left)[1] == (last1 - first1));
971  assert((M_bottom_right - M_upper_left)[0] == (result_last - result_first));
972  // typename RandomAccessIterator2d::difference_type size2d =
973 // M_bottom_right - M_upper_left;
974  typedef typename RandomAccessIterator2d::size_type size_type;
975  //size_type nb_rows = size2d[0];
976  //size_type nb_cols = size2d[1];
977  //typename std::iterator_traits<RandomAccessIterator1>::difference_type size1d = last1 - first1;
978  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
979 
980 
981 
982  size_type i = 0;
983  for(; result_first!=result_last; ++result_first, ++i)
984  {
985  *result_first = std::inner_product(M_upper_left.row_begin(i),
986  M_upper_left.row_end(i),
987  first1,
988  value_type());
989  }
990 
991  }
992 
993 
1026  template <class Matrix,
1027  class Vector1,
1028  class Vector2>
1029  inline
1030  void matrix_vector_multiplies(const Matrix& M,
1031  const Vector1& V,
1032  Vector2& Result)
1033  {
1034  slip::matrix_vector_multiplies(M.upper_left(),M.bottom_right(),
1035  V.begin(),V.end(),
1036  Result.begin(), Result.end());
1037  }
1038 
1074  template <class RandomAccessIterator2d,
1075  class RandomAccessIterator1,
1076  class RandomAccessIterator2>
1077  inline
1078  void hvector_matrix_multiplies(RandomAccessIterator1 V_first,
1079  RandomAccessIterator1 V_last,
1080  RandomAccessIterator2d M_up,
1081  RandomAccessIterator2d M_bot,
1082  RandomAccessIterator2 result_first,
1083  RandomAccessIterator2 result_last)
1084  {
1085  typename RandomAccessIterator2d::difference_type size2d =
1086  M_bot - M_up;
1087  typedef typename RandomAccessIterator2d::size_type size_type;
1088  size_type nb_rows = size2d[0];
1089  size_type nb_cols = size2d[1];
1090  typename std::iterator_traits<RandomAccessIterator1>::difference_type size1d =
1091  V_last - V_first;
1092 
1093  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
1094 
1095  assert(nb_rows == (std::size_t)size1d);
1096  assert(nb_cols == (std::size_t)(result_last - result_first));
1097 
1098  size_type j = 0;
1099  for(; result_first!=result_last; ++result_first, ++j)
1100  {
1101  *result_first = slip::hermitian_inner_product(V_first,
1102  V_last,
1103  M_up.col_begin(j),
1104  value_type());
1105  }
1106 
1107  }
1108 
1109 
1139  template <class Vector1,
1140  class Matrix,
1141  class Vector2>
1142  inline
1143  void hvector_matrix_multiplies(const Vector1& V,
1144  const Matrix& M,
1145  Vector2& R)
1146  {
1147  slip::hvector_matrix_multiplies(V.begin(),V.end(),
1148  M.upper_left(),M.bottom_right(),
1149  R.begin(),R.end());
1150  }
1185  template <class RandomAccessIterator2d,
1186  class RandomAccessIterator1,
1187  class RandomAccessIterator2>
1188  inline
1189  void tvector_matrix_multiplies(RandomAccessIterator1 V_first,
1190  RandomAccessIterator1 V_last,
1191  RandomAccessIterator2d M_up,
1192  RandomAccessIterator2d M_bot,
1193  RandomAccessIterator2 result_first,
1194  RandomAccessIterator2 result_last)
1195  {
1196 
1197  slip::hvector_matrix_multiplies(V_first,V_last,
1198  M_up,M_bot,
1199  result_first,result_last);
1200 
1201  }
1202 
1203 
1233  template <class Vector1,
1234  class Matrix,
1235  class Vector2>
1236  inline
1237  void tvector_matrix_multiplies(const Vector1& V,
1238  const Matrix& M,
1239  Vector2& R)
1240  {
1241  slip::tvector_matrix_multiplies(V.begin(),V.end(),
1242  M.upper_left(),M.bottom_right(),
1243  R.begin(),R.end());
1244  }
1245 
1246 
1247 
1285  template <class RandomAccessIterator2d1,
1286  class RandomAccessIterator2d2,
1287  class RandomAccessIterator2d3>
1288  inline
1289  void hmatrix_matrix_multiplies(RandomAccessIterator2d1 M1_up,
1290  RandomAccessIterator2d1 M1_bot,
1291  RandomAccessIterator2d2 M2_up,
1292  RandomAccessIterator2d2 M2_bot,
1293  RandomAccessIterator2d3 result_up,
1294  RandomAccessIterator2d3 result_bot)
1295  {
1296  assert((M1_bot-M1_up)[0] == (M2_bot-M2_up)[0]);
1297  assert((result_bot-result_up)[0] == (M1_bot-M1_up)[1]);
1298  assert((result_bot-result_up)[1] == (M2_bot-M2_up)[1]);
1299 
1300 
1301  typedef typename RandomAccessIterator2d1::size_type size_type;
1302  size_type rows = static_cast<size_type>((M1_bot-M1_up)[1]);
1303  size_type cols = static_cast<size_type>((M2_bot-M2_up)[1]);
1304 
1305 
1306  typedef typename std::iterator_traits<RandomAccessIterator2d3>::value_type value_type;
1307 
1308  slip::DPoint2d<int> d(0,0);
1309  for(size_type i = 0; i < rows; ++i)
1310  {
1311  d[0] = i;
1312  for(size_type j = 0; j < cols; ++j)
1313  {
1314  d[1] = j;
1315  result_up[d] =
1316  slip::hermitian_inner_product(M1_up.col_begin(i),M1_up.col_end(i),
1317  M2_up.col_begin(j),
1318  value_type());
1319  }
1320  }
1321  }
1322 
1355  template <class Matrix1,
1356  class Matrix2,
1357  class Matrix3>
1358  inline
1359  void hmatrix_matrix_multiplies(const Matrix1& M1,
1360  const Matrix2& M2,
1361  Matrix3& Result)
1362  {
1363  slip::hmatrix_matrix_multiplies(M1.upper_left(),M1.bottom_right(),
1364  M2.upper_left(),M2.bottom_right(),
1365  Result.upper_left(),Result.bottom_right());
1366  }
1367 
1368 
1399  template <class RandomAccessIterator2d,
1400  class RandomAccessIterator1,
1401  class RandomAccessIterator2>
1402  inline
1403  void hmatrix_vector_multiplies(RandomAccessIterator2d M_up,
1404  RandomAccessIterator2d M_bot,
1405  RandomAccessIterator1 V_first,
1406  RandomAccessIterator1 V_last,
1407  RandomAccessIterator2 R_first,
1408  RandomAccessIterator2 R_last)
1409  {
1410  assert((M_bot-M_up)[0] == (V_last-V_first));
1411 
1412  typedef typename RandomAccessIterator2d::size_type size_type;
1413  size_type cols = static_cast<size_type>((M_bot-M_up)[1]);
1414 
1415 
1416  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
1417 
1418  for(size_type j = 0; j < cols; ++j)
1419  {
1420 
1421  *R_first++ =
1422  slip::hermitian_inner_product(M_up.col_begin(j),
1423  M_up.col_end(j),
1424  V_first,
1425  value_type());
1426 
1427 
1428  }
1429  }
1430 
1431 
1456  template <class Matrix,
1457  class Vector1,
1458  class Vector2>
1459  inline
1460  void hmatrix_vector_multiplies(const Matrix& M,
1461  const Vector1& V,
1462  Vector2& R)
1463  {
1464  slip::hmatrix_vector_multiplies(M.upper_left(),M.bottom_right(),
1465  V.begin(),V.end(),
1466  R.begin(),R.end());
1467  }
1468 
1506  template <typename RandomAccessIterator2d,
1507  typename RandomAccessIterator1,
1508  typename RandomAccessIterator2>
1509  inline
1510  void gen_matrix_vector_multiplies(RandomAccessIterator2d M_up,
1511  RandomAccessIterator2d M_bot,
1512  RandomAccessIterator1 V_first,
1513  RandomAccessIterator1 V_last,
1514  RandomAccessIterator2 Y_first,
1515  RandomAccessIterator2 Y_last)
1516  {
1517  assert((M_bot - M_up)[1] == (V_last - V_first));
1518  assert((Y_last - Y_first) == (M_bot - M_up)[0]);
1519 
1520  typedef typename std::iterator_traits<RandomAccessIterator2d>::value_type value_type;
1521  typedef typename RandomAccessIterator2d::size_type size_type;
1522  typedef RandomAccessIterator2 iterator;
1523 
1524  iterator it = Y_first;
1525  iterator ite = Y_last;
1526  size_type i = 0;
1527  for(; it!=ite; ++it, ++i)
1528  {
1529  *it = *it + std::inner_product(M_up.row_begin(i),M_up.row_end(i),
1530  V_first,
1531  value_type());
1532  }
1533 
1534  }
1535 
1570  template <class Matrix,
1571  class Vector1,
1572  class Vector2>
1573  inline
1574  void gen_matrix_vector_multiplies(const Matrix& M,
1575  const Vector1& V,
1576  Vector2& Y)
1577  {
1578  slip::gen_matrix_vector_multiplies(M.upper_left(),M.bottom_right(),
1579  V.begin(),V.end(),
1580  Y.begin(),Y.end());
1581  }
1582 
1583 
1584 
1616  template <typename MatrixIterator1,
1617  typename MatrixIterator2,
1618  typename MatrixIterator3>
1619  inline
1620  void
1621  matrix_matrix_multiplies(MatrixIterator1 M1_up, MatrixIterator1 M1_bot,
1622  MatrixIterator2 M2_up, MatrixIterator2 M2_bot,
1623  MatrixIterator3 R_up, MatrixIterator3 R_bot)
1624 
1625  {
1626  assert((M1_bot - M1_up)[1] == (M2_bot - M2_up)[0]);
1627  assert((M1_bot - M1_up)[0] == (R_bot - R_up)[0]);
1628  assert((M2_bot - M2_up)[1] == (R_bot - R_up)[1]);
1629  typedef typename MatrixIterator1::value_type value_type;
1630  typedef typename MatrixIterator1::size_type size_type;
1631  size_type M1_rows = static_cast<size_type>((M1_bot - M1_up)[0]);
1632  size_type M2_cols = static_cast<size_type>((M2_bot - M2_up)[1]);
1633 
1634  for(size_type i = 0; i < M1_rows; ++i)
1635  {
1636  for(size_type j = 0; j < M2_cols; ++j)
1637  {
1638 
1639  *R_up = std::inner_product(M1_up.row_begin(i),M1_up.row_end(i),
1640  M2_up.col_begin(j),
1641  value_type());
1642  ++R_up;
1643  }
1644  }
1645 
1646  }
1647 
1648 
1684  template <class Matrix1,
1685  class Matrix2,
1686  class Matrix3>
1687  inline
1688  void matrix_matrix_multiplies(const Matrix1& M1,
1689  const Matrix2& M2,
1690  Matrix3& Result)
1691  {
1692  slip::matrix_matrix_multiplies(M1.upper_left(),M1.bottom_right(),
1693  M2.upper_left(),M2.bottom_right(),
1694  Result.upper_left(),Result.bottom_right());
1695  }
1696 
1697 
1737  template <typename RandomAccessIterator2d1,
1738  typename RandomAccessIterator2d2,
1739  typename RandomAccessIterator2d3>
1740  inline
1741  void gen_matrix_matrix_multiplies(RandomAccessIterator2d1 A_up,
1742  RandomAccessIterator2d1 A_bot,
1743  RandomAccessIterator2d2 B_up,
1744  RandomAccessIterator2d2 B_bot,
1745  RandomAccessIterator2d3 C_up,
1746  RandomAccessIterator2d3 C_bot)
1747  {
1748  assert((A_bot-A_up)[1] == (B_bot-B_up)[0]);
1749  assert((A_bot-A_up)[0] == (C_bot-C_up)[0]);
1750  assert((B_bot-B_up)[1] == (C_bot-C_up)[1]);
1751 
1752 
1753  typedef typename std::iterator_traits<RandomAccessIterator2d1>::value_type value_type;
1754  typedef typename RandomAccessIterator2d1::size_type size_type;
1755  size_type A_rows = (A_bot-A_up)[0];
1756  size_type B_cols = (B_bot-B_up)[1];
1757  for(size_type i = 0; i < A_rows; ++i)
1758  {
1759  for(size_type j = 0; j < B_cols; ++j)
1760  {
1761  *C_up += std::inner_product(A_up.row_begin(i),A_up.row_end(i),
1762  B_up.col_begin(j),
1763  value_type());
1764  ++C_up;
1765  }
1766  }
1767  }
1768 
1807  template <class Matrix1,
1808  class Matrix2,
1809  class Matrix3>
1810  inline
1811  void gen_matrix_matrix_multiplies(const Matrix1& A,
1812  const Matrix2& B,
1813  Matrix3& C)
1814  {
1815  slip::gen_matrix_matrix_multiplies(A.upper_left(),A.bottom_right(),
1816  B.upper_left(),B.bottom_right(),
1817  C.upper_left(),C.bottom_right());
1818  }
1819 
1820 
1844  template <class Matrix1, class Matrix2, class Matrix3>
1845  inline
1846  void multiplies(const Matrix1 & M1,
1847  const std::size_t nr1,
1848  const std::size_t nc1,
1849  const Matrix2 & M2,
1850  const std::size_t nr2,
1851  const std::size_t nc2,
1852  Matrix3 & Result,
1853  const std::size_t nr3,
1854  const std::size_t nc3)
1855  {
1856  assert(nc1 == nr2);
1857  assert(nr1 == nr3);
1858  assert(nc2 == nc3);
1859 
1860  for(std::size_t i = 0; i < nr1; ++i)
1861  {
1862  for(std::size_t j = 0; j < nc2; ++j)
1863  {
1864  Result[i][j] = 0;
1865  for(std::size_t k = 0; k < nc1; ++k)
1866  {
1867  Result[i][j] += M1[i][k] * M2[k][j];
1868  }
1869  }
1870  }
1871  }
1872 
1894  template <class Matrix1, class Matrix2, class Matrix3>
1895  inline
1896  void multiplies(const Matrix1 & M,
1897  const std::size_t nrm,
1898  const std::size_t ncm,
1899  const Matrix2 & V,
1900  const std::size_t nrv,
1901  Matrix3 & Result,
1902  const std::size_t nrr)
1903  {
1904  assert(ncm == nrv);
1905  assert(nrm == nrr);
1906 
1907  for(std::size_t i = 0; i < nrm; ++i)
1908  {
1909  Result[i] = 0;
1910  for(std::size_t k = 0; k < ncm; ++k)
1911  {
1912  Result[i] += M[i][k] * V[k];
1913  }
1914  }
1915  }
1916 
1917 
1948  template <class Matrix1, class Matrix2, class Matrix3>
1949  inline
1950  void multiplies(const Matrix1 & M1,
1951  const Matrix2 & M2,
1952  Matrix3 & Result)
1953  {
1954  assert(M1.cols() == M2.rows());
1955  assert(M1.rows() == Result.rows());
1956  assert(M2.cols() == Result.cols());
1957  std::size_t nr1 = M1.rows();
1958  std::size_t nc1 = M1.cols();
1959  std::size_t nc2 = M2.cols();
1960  slip::multiplies(M1,nr1,nc1,M2,nc1,nc2,Result,nr1,nc2);
1961  }
1962 
1983  template <typename MatrixIterator1,
1984  typename MatrixIterator2,
1985  typename MatrixIterator3>
1986  inline
1987  void multiplies(MatrixIterator1 A1_up, MatrixIterator1 A1_bot,
1988  MatrixIterator2 A2_up, MatrixIterator2 A2_bot,
1989  MatrixIterator3 Res_up, MatrixIterator3 Res_bot)
1990  {
1991  typedef typename std::iterator_traits<MatrixIterator1>::value_type _T;
1992  typename MatrixIterator1::difference_type sizeA1= A1_bot - A1_up;
1993  typename MatrixIterator2::difference_type sizeA2= A2_bot - A2_up;
1994  typename MatrixIterator3::difference_type sizeRes= Res_bot - Res_up;
1995  assert(sizeRes[0] == sizeA1[0]);
1996  assert(sizeRes[1] == sizeA2[1]);
1997  assert(sizeA1[1] == sizeA2[0]);
1998 
1999  for(int i = 0; i < sizeRes[0]; ++i)
2000  {
2001  for(int j = 0; j < sizeRes[1]; ++j)
2002  {
2003  slip::DPoint2d<int> ij(i,j);
2004  *(Res_up + ij) = _T(0);
2005  for(int k = 0; k < sizeA1[1]; ++k)
2006  {
2007  slip::DPoint2d<int> ik(i,k);
2008  slip::DPoint2d<int> kj(k,j);
2009  *(Res_up + ij) += (*(A1_up + ik)) * (*(A2_up + kj));
2010  }
2011  }
2012  }
2013  }
2014 
2044  template <class RandomAccessIterator2d1,
2045  class RandomAccessIterator2d2>
2046  inline
2047  void matrix_matrixt_multiplies(RandomAccessIterator2d1 M1_up,
2048  RandomAccessIterator2d1 M1_bot,
2049  RandomAccessIterator2d2 Res_up,
2050  RandomAccessIterator2d2 Res_bot)
2051  {
2052  assert((Res_bot-Res_up)[0] == (M1_bot-M1_up)[0]);
2053  assert((Res_bot-Res_up)[1] == (M1_bot-M1_up)[0]);
2054 
2055 
2056  typedef typename RandomAccessIterator2d1::size_type size_type;
2057  size_type rows = static_cast<size_type>((M1_bot-M1_up)[0]);
2058 
2059  typedef typename std::iterator_traits<RandomAccessIterator2d2>::value_type value_type;
2060 
2061  slip::DPoint2d<int> d(0,0);
2062  for(size_type i = 0; i < rows; ++i)
2063  {
2064  d[0] = i;
2065  for(size_type j = 0; j < rows; ++j)
2066  {
2067  d[1] = j;
2068  Res_up[d] =
2069  std::inner_product(M1_up.row_begin(i),M1_up.row_end(i),
2070  M1_up.row_begin(j),
2071  value_type());
2072  }
2073  }
2074  }
2075 
2076 
2100  template <class Matrix1,
2101  class Matrix2>
2102  inline
2103  void matrix_matrixt_multiplies(const Matrix1& M1,
2104  Matrix2& result)
2105  {
2106  slip::matrix_matrixt_multiplies(M1.upper_left(),M1.bottom_right(),
2107  result.upper_left(),result.bottom_right());
2108  }
2109 
2142  template <class RandomAccessIterator2d1,
2143  class T,
2144  class RandomAccessIterator2d2>
2145  inline
2146  void matrix_shift(RandomAccessIterator2d1 M_upper_left,
2147  RandomAccessIterator2d1 M_bottom_right,
2148  const T& lambda,
2149  RandomAccessIterator2d2 R_upper_left,
2150  RandomAccessIterator2d2 R_bottom_right)
2151  {
2152  assert((M_bottom_right - M_upper_left)[0] == (M_bottom_right - M_upper_left)[1]);
2153  assert((M_bottom_right - M_upper_left)[0] == (R_bottom_right - R_upper_left)[0]);
2154  assert((M_bottom_right - M_upper_left)[1] == (R_bottom_right - R_upper_left)[1]);
2155  typename RandomAccessIterator2d1::difference_type size2d_M =
2156  M_bottom_right - M_upper_left;
2157  typename RandomAccessIterator2d2::difference_type size2d_R =
2158  R_bottom_right - R_upper_left;
2159  typedef typename RandomAccessIterator2d1::size_type size_type;
2160  size_type M_rows = size2d_M[0];
2161 
2162 
2163 
2164  typedef typename std::iterator_traits<RandomAccessIterator2d2>::value_type value_type;
2165 
2166  std::copy(M_upper_left,M_bottom_right,R_upper_left);
2167  slip::DPoint2d<int> d(1,1);
2168  for(size_type i = 0; i < M_rows; ++i)
2169  {
2170  *R_upper_left = value_type(*M_upper_left) - value_type(lambda);
2171  R_upper_left += d;
2172  M_upper_left += d;
2173  }
2174 
2175  }
2176 
2203  template <class RandomAccessIterator2d1,
2204  class T>
2205  inline
2206  void matrix_shift(RandomAccessIterator2d1 M_upper_left,
2207  RandomAccessIterator2d1 M_bottom_right,
2208  const T& lambda)
2209  {
2210  typename RandomAccessIterator2d1::difference_type size2d_M =
2211  M_bottom_right - M_upper_left;
2212  typedef typename RandomAccessIterator2d1::size_type size_type;
2213  size_type M_rows = size2d_M[0];
2214  size_type M_cols = size2d_M[1];
2215  assert(M_rows == M_cols);
2216 
2217  typedef typename std::iterator_traits<RandomAccessIterator2d1>::value_type value_type;
2218 
2219  slip::DPoint2d<int> d(1,1);
2220  for(size_type i = 0; i < M_rows; ++i)
2221  {
2222  *M_upper_left -= value_type(lambda);
2223  M_upper_left += d;
2224  }
2225 
2226  }
2227 
2255  template <class Matrix1,
2256  class T,
2257  class Matrix2>
2258  inline
2259  void matrix_shift(const Matrix1& A,
2260  const T& lambda,
2261  Matrix2& R)
2262  {
2263  slip::matrix_shift(A.upper_left(),A.bottom_right(),
2264  lambda,
2265  R.upper_left(),R.bottom_right());
2266 
2267  }
2268 
2292  template <class Matrix,
2293  class T>
2294  inline
2295  void matrix_shift(Matrix& A,
2296  const T& lambda)
2297  {
2298  slip::matrix_shift(A.upper_left(),A.bottom_right(),lambda);
2299 
2300  }
2301 
2302 
2341  template <typename RandomAccessIterator2d1,
2342  typename T,
2343  typename RandomAccessIterator1,
2344  typename RandomAccessIterator2>
2345  inline
2346  void rank1_update(RandomAccessIterator2d1 A_up,
2347  RandomAccessIterator2d1 A_bot,
2348  const T& alpha,
2349  RandomAccessIterator1 X_first,
2350  RandomAccessIterator1 X_last,
2351  RandomAccessIterator2 Y_first,
2352  RandomAccessIterator2 Y_last)
2353  {
2354 
2355 
2356  assert((A_bot - A_up)[0] == (X_last - X_first));
2357  assert((A_bot - A_up)[1] == (Y_last - Y_first));
2358 
2359  typename RandomAccessIterator2d1::difference_type size2d_A =
2360  A_bot - A_up;
2361  typedef typename RandomAccessIterator2d1::size_type size_type;
2362  size_type A_rows = size2d_A[0];
2363 
2364  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference;
2365  _Difference y_size = Y_last - Y_first;
2366 
2367 
2368  typedef typename std::iterator_traits<RandomAccessIterator2d1>::value_type value_type;
2369 
2370 
2371 
2372  slip::Array<value_type> tmp((std::size_t)(y_size));
2373  //computes tmp = Y^HA
2374  for(size_type i = 0; i < A_rows; ++i, ++X_first)
2375  {
2376  slip::conj_vector_scalar_multiplies(Y_first,Y_last,
2377  *X_first,
2378  tmp.begin());
2379  slip::aXpY(alpha,tmp.begin(),tmp.end(),A_up.row_begin(i));
2380  }
2381  }
2382 
2383 
2416  template <typename Matrix,
2417  typename T,
2418  typename Vector1,
2419  typename Vector2>
2420  inline
2421  void rank1_update(Matrix& A,
2422  const T& alpha,
2423  const Vector1& X,
2424  const Vector2& Y)
2425  {
2426  slip::rank1_update(A.upper_left(),A.bottom_right(),
2427  alpha,
2428  X.begin(),X.end(),
2429  Y.begin(),Y.end());
2430  }
2431 
2446  template <typename Iterator1, typename Iterator2, typename U>
2447  inline
2448  void reduction(Iterator1 seq1_beg, Iterator1 seq1_end,
2449  Iterator2 seq2_beg, U u)
2450  {
2451  typedef typename std::iterator_traits<Iterator1>::value_type _T1;
2452  typedef typename std::iterator_traits<Iterator2>::value_type _T2;
2453  std::transform(seq2_beg,seq2_beg + (seq1_end - seq1_beg),
2454  seq1_beg,
2455  seq1_beg,
2456  slip::saxpy<_T2,_T1,_T1,U>(u));
2457  }
2458 
2474  template <class Matrix1, class Matrix2>
2475  inline
2476  void hermitian_transpose(const Matrix1 & M, Matrix2 & TM)
2477  {
2478  assert(M.rows() == TM.cols());
2479  assert(M.cols() == TM.rows());
2480  typedef typename Matrix1::value_type _T;
2481  std::size_t M_rows = M.rows();
2482  std::size_t M_cols = M.cols();
2483  for (std::size_t i = 0; i < M_rows; ++i)
2484  {
2485  for (std::size_t j = 0; j < M_cols; ++j)
2486  {
2487  TM[j][i] = slip::conj(M[i][j]);
2488  }
2489  }
2490 
2491  }
2492 
2493 
2494 
2514  template <class Matrix1, class Matrix2>
2515  inline
2516  void transpose(const Matrix1 & M,
2517  const std::size_t nr1,
2518  const std::size_t nc1,
2519  Matrix2 & TM,
2520  const std::size_t nr2,
2521  const std::size_t nc2)
2522  {
2523  assert(nr1==nc2);
2524  assert(nc1==nr2);
2525 
2526  for (std::size_t i = 0; i < nr1; ++i)
2527  {
2528  for (std::size_t j = 0; j < nc1; ++j)
2529  {
2530  TM[j][i] = M[i][j];
2531  }
2532  }
2533 
2534  }
2535 
2560  template <class Matrix1, class Matrix2>
2561  inline
2562  void transpose(const Matrix1 & M, Matrix2 & TM)
2563  {
2564 
2565  slip::transpose(M,M.rows(),M.cols(),TM,TM.rows(),TM.cols());
2566  }
2567 
2568 
2601  template <class Matrix1, class Matrix2>
2602  inline
2603  void conj(const Matrix1 & M, Matrix2 & CM)
2604  {
2605  assert(M.rows() == CM.rows());
2606  assert(M.cols() == CM.cols());
2607  slip::apply(M.begin(),M.end(),CM.begin(),slip::conj);
2608  }
2609 
2641  template <class Matrix1, class Matrix2>
2642  inline
2643  void real(const Matrix1 & C, Matrix2 & R)
2644  {
2645  assert(R.rows() == C.rows());
2646  assert(R.cols() == C.cols());
2647  slip::apply(C.begin(),C.end(),R.begin(),std::real);
2648  }
2649 
2650 
2682  template <class Matrix1, class Matrix2>
2683  inline
2684  void imag(const Matrix1 & C, Matrix2 & I)
2685  {
2686  assert(I.rows() == C.rows());
2687  assert(I.cols() == C.cols());
2688  slip::apply(C.begin(),C.end(),I.begin(),std::imag);
2689  }
2690 
2691 
2692 
2693 
2694 
2735  template <typename VectorIterator1,
2736  typename VectorIterator2,
2737  typename MatrixIterator>
2738  inline
2739  void outer_product(VectorIterator1 V_begin, VectorIterator1 V_end,
2740  VectorIterator2 W_begin, VectorIterator2 W_end,
2741  MatrixIterator R_up, MatrixIterator R_bot)
2742  {
2743 
2744  typename MatrixIterator::difference_type sizeR = R_bot - R_up;
2745  std::size_t R_rows = sizeR[0];
2746 
2747  assert(sizeR[0] == (V_end - V_begin));
2748  assert(sizeR[1] == (W_end - W_begin));
2749 
2750  for(std::size_t i = 0; i < R_rows; ++i, ++V_begin)
2751  {
2752  slip::conj_vector_scalar_multiplies(W_begin,W_end,
2753  *V_begin,
2754  R_up.row_begin(i));
2755  }
2756  }
2757 
2856 template<class RandomAccessIterator1,
2857  class RandomAccessIterator2,
2858  class RandomAccessIterator2d>
2859  void rank1_tensorial_product(RandomAccessIterator1 base1_first,
2860  RandomAccessIterator1 base1_end,
2861  RandomAccessIterator2 base2_first,
2862  RandomAccessIterator2 base2_end,
2863  RandomAccessIterator2d matrix_upper_left,
2864  RandomAccessIterator2d matrix_bottom_right)
2865 {
2866  assert((matrix_bottom_right - matrix_upper_left)[0] == (base1_end - base1_first));
2867  assert((matrix_bottom_right - matrix_upper_left)[1] == (base2_end - base2_first));
2868 
2869  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Distance1;
2870  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance2;
2871 
2872  typedef typename std::iterator_traits<RandomAccessIterator2d>::difference_type _Distance2d;
2873 
2874  _Distance1 base1_first_size = (base1_end - base1_first);
2875  _Distance2 base2_first_size = (base2_end - base2_first);
2876 
2877  _Distance2d matrix_size2d = (matrix_bottom_right - matrix_upper_left);
2878 
2879  for(_Distance2 i = 0 ; i < base1_first_size; ++i)
2880  {
2881  for(_Distance1 j = 0 ; j < base2_first_size; ++j)
2882  {
2883  *(matrix_upper_left++) = base1_first[i] * base2_first[j];
2884  }
2885  }
2886 
2887 }
2888 
2889 
3033 template<class RandomAccessIterator1,
3034  class RandomAccessIterator2,
3035  class RandomAccessIterator3,
3036  class RandomAccessIterator3d>
3037  void rank1_tensorial_product(RandomAccessIterator1 base1_first,
3038  RandomAccessIterator1 base1_end,
3039  RandomAccessIterator2 base2_first,
3040  RandomAccessIterator2 base2_end,
3041  RandomAccessIterator3 base3_first,
3042  RandomAccessIterator3 base3_end,
3043  RandomAccessIterator3d base_front_upper_left,
3044  RandomAccessIterator3d base_back_bottom_right)
3045 {
3046  assert((base_back_bottom_right - base_front_upper_left)[0] == (base1_end - base1_first));
3047  assert((base_back_bottom_right - base_front_upper_left)[1] == (base2_end - base2_first));
3048  assert((base_back_bottom_right - base_front_upper_left)[2] == (base3_end - base3_first));
3049 
3050  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Distance1;
3051  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance2;
3052  typedef typename std::iterator_traits<RandomAccessIterator3>::difference_type _Distance3;
3053 
3054  typedef typename std::iterator_traits<RandomAccessIterator3d>::difference_type _Distance3d;
3055 
3056  _Distance1 base1_first_size = (base1_end - base1_first);
3057  _Distance2 base2_first_size = (base2_end - base2_first);
3058  _Distance2 base3_first_size = (base3_end - base3_first);
3059 
3060  _Distance3d base_size3d = (base_back_bottom_right - base_front_upper_left);
3061  for(_Distance3 k = 0 ; k < base1_first_size; ++k)
3062  {
3063  for(_Distance2 i = 0 ; i < base2_first_size; ++i)
3064  {
3065  for(_Distance1 j = 0 ; j < base3_first_size; ++j)
3066  {
3067  *(base_front_upper_left++) = base1_first[k]*base2_first[i] * base3_first[j];
3068  }
3069  }
3070  }
3071 
3072 }
3073 
3074 
3075  /* @} */
3076 
3077 
3078 
3080  /* @{ */
3081 
3102  template<typename InputIterator>
3103  inline
3104  typename slip::lin_alg_traits<typename std::iterator_traits<InputIterator>::value_type>::value_type
3105  L22_norm_cplx(InputIterator first,
3106  InputIterator last)
3107 
3108 {
3109  typedef typename std::iterator_traits<InputIterator>::value_type _T;
3110  _T __init = _T();
3111  for (; first != last; ++first)
3112  {
3113  __init +=slip::conj(*first)*(*first);
3114  }
3115  return std::real(__init);
3116 }
3117 
3118 
3139  template <typename RandomAccessIterator2d>
3140  inline
3141  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator2d>::value_type>::value_type row_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
3142  {
3143  assert(((bottom_right - upper_left)[0] != 0) && ((bottom_right - upper_left)[1] != 0));
3144 
3145  typedef typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator2d>::value_type>::value_type value_type;
3146  typedef typename RandomAccessIterator2d::size_type size_type;
3147  size_type rows = (bottom_right - upper_left)[0];
3148  value_type max = slip::L1_norm<value_type>(upper_left.row_begin(0),
3149  upper_left.row_end(0));
3150 
3151  for(std::size_t i = 1; i < rows; ++i)
3152  {
3153  value_type tmp = slip::L1_norm<value_type>(upper_left.row_begin(i),
3154  upper_left.row_end(i));
3155  if(tmp > max)
3156  {
3157  max = tmp;
3158  }
3159  }
3160  return max;
3161  }
3162 
3163 
3164 
3184  template <typename Container2d>
3185  inline
3186  typename slip::lin_alg_traits<typename Container2d::value_type>::value_type row_norm(const Container2d& container)
3187  {
3188  return slip::row_norm(container.upper_left(),container.bottom_right());
3189  }
3190 
3191 
3212  template <typename RandomAccessIterator2d>
3213  inline
3214  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator2d>::value_type>::value_type col_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
3215  {
3216  assert(((bottom_right - upper_left)[0] != 0) && ((bottom_right - upper_left)[1] != 0));
3217 
3218  typedef typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator2d>::value_type>::value_type value_type;
3219  typedef typename RandomAccessIterator2d::size_type size_type;
3220  size_type cols = (bottom_right - upper_left)[1];
3221  value_type max = slip::L1_norm<value_type>(upper_left.col_begin(0),
3222  upper_left.col_end(0));
3223 
3224  for(std::size_t j = 1; j < cols; ++j)
3225  {
3226  value_type tmp = slip::L1_norm<value_type>(upper_left.col_begin(j),
3227  upper_left.col_end(j));
3228  if(tmp > max)
3229  {
3230  max = tmp;
3231  }
3232  }
3233  return max;
3234  }
3235 
3236 
3256  template <typename Container2d>
3257  inline
3258  typename slip::lin_alg_traits<typename Container2d::value_type>::value_type col_norm(const Container2d& container)
3259  {
3260  return slip::col_norm(container.upper_left(),container.bottom_right());
3261  }
3262 
3263 
3284  template <typename RandomAccessIterator2d>
3285  inline
3286  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator2d>::value_type>::value_type frobenius_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
3287  {
3288  assert(((bottom_right - upper_left)[0] != 0) && ((bottom_right - upper_left)[1] != 0));
3289 
3290  return std::sqrt(slip::L22_norm_cplx(upper_left,bottom_right));
3291  }
3292 
3313  template <typename Container2d>
3314  inline
3315  typename slip::lin_alg_traits<typename Container2d::value_type>::value_type frobenius_norm(const Container2d& container)
3316  {
3317  return slip::frobenius_norm(container.upper_left(),container.bottom_right());
3318  }
3319 
3320  /* @} */
3321 
3322 
3323 
3324 
3325 
3327  /* @{ */
3328 
3368  template<typename Container2d,
3369  typename RandomAccessIterator1>
3370  void get_diagonal(const Container2d& container,
3371  RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
3372  const int diag_number = 0)
3373 
3374  {
3375  assert(container.dim1() == container.dim2());
3376  assert(std::abs(diag_number) < int(container.dim1()));
3377  assert((diag_last - diag_first) <= int(container.dim1()));
3378  int i = 0;
3379  int j = 0;
3380  if(diag_number >= 0)
3381  {
3382  i = 0;
3383  j = diag_number;
3384  }
3385  else
3386  {
3387  i = -diag_number;
3388  j = 0;
3389  }
3390 
3391  for(; diag_first != diag_last; ++diag_first)
3392  {
3393  *diag_first = container[i][j];
3394  i = i + 1;
3395  j = j + 1;
3396  }
3397  }
3398 
3440  template<typename Container2d,
3441  typename RandomAccessIterator1>
3442  void set_diagonal(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
3443  Container2d& container,
3444  const int diag_number = 0)
3445 
3446  {
3447  assert(container.dim1() == container.dim2());
3448  assert(std::abs(diag_number) < int(container.dim1()));
3449  assert((diag_last - diag_first) <= int(container.dim1()));
3450  int i = 0;
3451  int j = 0;
3452  if(diag_number >= 0)
3453  {
3454  i = 0;
3455  j = diag_number;
3456  }
3457  else
3458  {
3459  i = -diag_number;
3460  j = 0;
3461  }
3462 
3463  for(; diag_first != diag_last; ++diag_first)
3464  {
3465  container[i][j] = *diag_first;
3466  i = i + 1;
3467  j = j + 1;
3468  }
3469  }
3470 
3502  template<typename Container2d>
3503  void fill_diagonal(Container2d& container,
3504  const typename Container2d::value_type& val,
3505  const int diag_number = 0)
3506 
3507  {
3508  assert(container.dim1() == container.dim2());
3509  assert(std::abs(diag_number) < int(container.dim1()));
3510 
3511  int i = 0;
3512  int j = 0;
3513  if(diag_number >= 0)
3514  {
3515  i = 0;
3516  j = diag_number;
3517  }
3518  else
3519  {
3520  i = -diag_number;
3521  j = 0;
3522  }
3523 
3524  int n = int(container.dim1()) - std::abs(diag_number);
3525  for(int k = 0; k < n; ++k)
3526  {
3527  container[i][j] = val;
3528  i = i + 1;
3529  j = j + 1;
3530  }
3531  }
3532 
3555  template<typename MatrixIterator>
3556  void fill_diagonal(MatrixIterator M_up,
3557  MatrixIterator M_bot,
3558  const typename MatrixIterator::value_type& val,
3559  const int diag_number = 0)
3560 
3561  {
3562  //assert((M_bot-M_up)[0] == (M_bot-M_up)[1]);
3563  assert(std::abs(diag_number) < int((M_bot-M_up)[0]));
3564 
3565  int i = 0;
3566  int j = 0;
3567  if(diag_number >= 0)
3568  {
3569  i = 0;
3570  j = diag_number;
3571  }
3572  else
3573  {
3574  i = -diag_number;
3575  j = 0;
3576  }
3577 
3578  int n = int((M_bot-M_up)[0]) - std::abs(diag_number);
3579  slip::DPoint2d<int> di(i,j);
3580  slip::DPoint2d<int> d(1,1);
3581  for(int k = 0; k < n; ++k)
3582  {
3583  M_up[di] = val;
3584  di+= d;
3585  }
3586  }
3587 
3588 /* @} */
3589 
3591 /* @{ */
3606  template <typename Matrix>
3607  inline
3608  Matrix identity(const std::size_t nr, const std::size_t nc)
3609  {
3610  typedef typename Matrix::value_type _T;
3611  Matrix Res(nr,nc,_T(0));
3612  for(std::size_t i = 0; i < nr; ++i)
3613  {
3614  if(i < nc)
3615  {
3616  Res[i][i] = _T(1);
3617  }
3618  }
3619  return Res;
3620  }
3621 
3637  template <typename MatrixIterator>
3638  inline
3639  void identity(MatrixIterator A_up,
3640  MatrixIterator A_bot)
3641  {
3642  typedef typename MatrixIterator::value_type value_type;
3643  int A_rows = int((A_bot - A_up)[0]);
3644  int A_cols = int((A_bot - A_up)[1]);
3645  int n = std::min(A_rows,A_cols);
3646  std::fill(A_up,A_bot,value_type(0));
3647  slip::DPoint2d<int> d(1,1);
3648  for (int i = 0; i < n; ++i)
3649  {
3650  *A_up = value_type(1);
3651  A_up += d;
3652  }
3653  }
3654 
3670  template <typename Matrix>
3671  inline
3672  void identity(Matrix & A)
3673  {
3674  typedef typename Matrix::value_type _T;
3675  std::size_t A_rows = A.rows();
3676  std::size_t A_cols = A.cols();
3677  for(std::size_t i = 0; i < A_rows; ++i)
3678  {
3679  if(i < A_cols)
3680  {
3681  for(std::size_t j = 0; j < i; ++j)
3682  {
3683  A[i][j] = _T(0);
3684  }
3685  A[i][i] = _T(1);
3686  for(std::size_t j = i+1; j < A_cols; ++j)
3687  {
3688  A[i][j] = _T(0);
3689  }
3690  }
3691  else
3692  {
3693  for(std::size_t j = 0; j < A_cols; ++j)
3694  {
3695  A[i][j] = _T(0);
3696  }
3697  }
3698  }
3699  }
3700 
3715  template <typename Matrix>
3716  inline
3717  void hilbert(Matrix & A)
3718  {
3719  typedef typename Matrix::value_type value_type;
3720  std::size_t A_rows = A.rows();
3721  std::size_t A_cols = A.cols();
3722  for(std::size_t i = 0; i < A_rows; ++i)
3723  {
3724  for(std::size_t j = 0; j < A_cols; ++j)
3725  {
3726  A(i,j) = value_type(1) / value_type((i + j) + 1);
3727  }
3728  }
3729 
3730  }
3731 
3756  template <typename RandomAccessIterator, typename Matrix>
3757  inline
3758  void toeplitz(RandomAccessIterator first, RandomAccessIterator last,
3759  Matrix& A)
3760  {
3761  typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
3762  typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Difference;
3763 
3764  _Difference size = last -first;
3765  int n = int(size);
3766 
3767  assert(A.rows() == A.cols());
3768  assert(int(size) == int(A.rows()));
3769 
3770  slip::fill_diagonal(A,*first++);
3771 
3772  for(int i = 1; i < n; ++i, ++first)
3773  {
3774  slip::fill_diagonal(A,*first,i);
3775  slip::fill_diagonal(A,*first,-i);
3776  }
3777  }
3778 
3779 /* @} */
3780 
3782  /* @{ */
3798  template <typename Matrix>
3799  inline
3800  bool is_squared(const Matrix & A)
3801  {
3802  return (A.dim1() == A.dim2());
3803  }
3804 
3821  template <typename MatrixIterator1>
3822  inline
3823  bool is_squared(MatrixIterator1 A_up,
3824  MatrixIterator1 A_bot)
3825  {
3826  return ((A_bot - A_up)[0] == (A_bot - A_up)[1]);
3827  }
3828 
3846  template <typename Matrix>
3847  inline
3848  bool is_identity(const Matrix & A,
3849  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
3850  = slip::epsilon<typename Matrix::value_type>())
3851  {
3852  if(!is_squared(A))
3853  {
3854  return false;
3855  }
3856 
3857  typedef typename Matrix::value_type value_type;
3858  for(std::size_t i = 0; i < A.rows(); ++i)
3859  {
3860  if(std::abs(A[i][i] - value_type(1)) > tol)
3861  {
3862  return false;
3863  }
3864  for(std::size_t j = 0; j < i; ++j)
3865  {
3866  if((std::abs(A[i][j]) > tol) || (std::abs(A[j][i]) > tol))
3867  {
3868  return false;
3869  }
3870  }
3871  }
3872  return true;
3873 
3874  }
3875 
3876 
3877 
3897  template <typename MatrixIterator1>
3898  inline
3899  bool is_diagonal(MatrixIterator1 A_up,
3900  MatrixIterator1 A_bot,
3901  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
3902  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
3903  {
3904  if(!is_squared(A_up,A_bot))
3905  {
3906  return false;
3907  }
3908 
3909  typedef typename std::iterator_traits<MatrixIterator1>::value_type value_type;
3910  const std::size_t A_rows = (std::size_t)(A_bot - A_up)[0];
3911  slip::DPoint2d<int> d1;
3912  slip::DPoint2d<int> d2;
3913 
3914  for(std::size_t i = 0; i < A_rows; ++i)
3915  {
3916  for(std::size_t j = 0; j < i; ++j)
3917  {
3918  d1.set_coord(i,j);
3919  d2.set_coord(j,i);
3920  if((std::abs(A_up[d1]) > tol) || (std::abs(A_up[d2]) > tol))
3921  {
3922  return false;
3923  }
3924  }
3925  }
3926  return true;
3927 
3928  }
3929 
3930 
3931 
3950  template <typename Matrix>
3951  inline
3952  bool is_diagonal(const Matrix & A,
3953  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
3954  = slip::epsilon<typename Matrix::value_type>())
3955  {
3956 
3957  return slip::is_diagonal(A.upper_left(),A.bottom_right(),tol);
3958 
3959  }
3960 
3961 
3984  template <typename MatrixIterator1>
3985  inline
3986  bool is_hermitian(MatrixIterator1 A_up,
3987  MatrixIterator1 A_bot,
3988  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
3989  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
3990  {
3991  if(!is_squared(A_up,A_bot))
3992  {
3993  return false;
3994  }
3995 
3996  typedef typename std::iterator_traits<MatrixIterator1>::value_type value_type;
3997  std::size_t A_rows = (std::size_t)(A_bot - A_up)[0];
3998  slip::DPoint2d<int> d1;
3999  slip::DPoint2d<int> d2;
4000 
4001  for(std::size_t i = 0; i < A_rows; ++i)
4002  {
4003  for(std::size_t j = 0; j <= i; ++j)
4004  {
4005  d1.set_coord(i,j);
4006  d2.set_coord(j,i);
4007  if(std::abs(A_up[d1] - slip::conj(A_up[d2])) > tol)
4008  {
4009  return false;
4010  }
4011  }
4012  }
4013  return true;
4014 
4015  }
4016 
4038  template <typename Matrix>
4039  inline
4040  bool is_hermitian(const Matrix & A,
4041  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4042  = slip::epsilon<typename Matrix::value_type>())
4043  {
4044 
4045  return slip::is_hermitian(A.upper_left(),A.bottom_right(),tol);
4046 
4047  }
4082  template <typename MatrixIterator1>
4083  inline
4084  bool is_skew_hermitian(MatrixIterator1 A_up,
4085  MatrixIterator1 A_bot,
4086  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4087  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4088  {
4089  if(!is_squared(A_up,A_bot))
4090  {
4091  return false;
4092  }
4093 
4094  typedef typename std::iterator_traits<MatrixIterator1>::value_type value_type;
4095  std::size_t A_rows = (std::size_t)(A_bot - A_up)[0];
4096  slip::DPoint2d<int> d1;
4097  slip::DPoint2d<int> d2;
4098 
4099  for(std::size_t i = 0; i < A_rows; ++i)
4100  {
4101  for(std::size_t j = 0; j <= i; ++j)
4102  {
4103  d1.set_coord(i,j);
4104  d2.set_coord(j,i);
4105  if(std::abs(A_up[d1] + slip::conj(A_up[d2])) > tol)
4106  {
4107  return false;
4108  }
4109  }
4110  }
4111  return true;
4112 
4113  }
4114 
4148  template <typename Matrix>
4149  inline
4150  bool is_skew_hermitian(const Matrix & A,
4151  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4152  = slip::epsilon<typename Matrix::value_type>())
4153  {
4154 
4155  return slip::is_hermitian(A.upper_left(),A.bottom_right(),tol);
4156 
4157  }
4158 
4159 
4187  template <typename MatrixIterator1>
4188  inline
4189  bool is_symmetric(MatrixIterator1 A_up,
4190  MatrixIterator1 A_bot,
4191  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4192  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4193  {
4194  if(slip::is_complex(typename std::iterator_traits<MatrixIterator1>::value_type()))
4195  {
4196  return false;
4197  }
4198  return slip::is_hermitian(A_up,A_bot,tol);
4199  }
4200 
4227  template <typename Matrix>
4228  inline
4229  bool is_symmetric(const Matrix & A,
4230  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4231  = slip::epsilon<typename Matrix::value_type>())
4232  {
4233 
4234  return slip::is_symmetric(A.upper_left(),A.bottom_right(),tol);
4235 
4236  }
4237 
4238 
4239 
4240 
4264  template <typename MatrixIterator1>
4265  inline
4266  bool is_band_matrix(MatrixIterator1 A_up,
4267  MatrixIterator1 A_bot,
4268  const typename MatrixIterator1::size_type lower_band_width,
4269  const typename MatrixIterator1::size_type upper_band_width,
4270  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4271  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4272  {
4273 
4274  assert(typename MatrixIterator1::size_type(lower_band_width) < typename MatrixIterator1::size_type((A_bot - A_up)[0]));
4275  assert(typename MatrixIterator1::size_type(upper_band_width) < typename MatrixIterator1::size_type((A_bot - A_up)[0]));
4276  assert(lower_band_width >= typename MatrixIterator1::size_type(0));
4277  assert(upper_band_width >= typename MatrixIterator1::size_type(0));
4278 
4279  if(!is_squared(A_up,A_bot))
4280  {
4281  return false;
4282  }
4283  typedef typename std::iterator_traits<MatrixIterator1>::value_type value_type;
4284  std::size_t A_rows = (std::size_t)(A_bot - A_up)[0];
4285 
4286  slip::DPoint2d<int> d1;
4287 
4288  for(std::size_t i = 0; i < (A_rows - upper_band_width); ++i)
4289  {
4290  for(std::size_t j = (i + upper_band_width + 1); j < A_rows; ++j)
4291  {
4292  d1.set_coord(i,j);
4293  if(std::abs(A_up[d1]) > tol)
4294  {
4295  return false;
4296  }
4297  }
4298  }
4299 
4300  for(std::size_t i = (lower_band_width + 1); i < A_rows; ++i)
4301  {
4302  for(std::size_t j = 0; j < (i-lower_band_width); ++j)
4303  {
4304  d1.set_coord(i,j);
4305  if(std::abs(A_up[d1]) > tol)
4306  {
4307  return false;
4308  }
4309  }
4310  }
4311  return true;
4312  }
4313 
4331  template <typename Matrix>
4332  inline
4333  bool is_band_matrix(const Matrix & A,
4334  const typename Matrix::size_type lower_band_width,
4335  const typename Matrix::size_type upper_band_width,
4336  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4337  = slip::epsilon<typename Matrix::value_type>())
4338  {
4339 
4340  return slip::is_band_matrix(A.upper_left(),A.bottom_right(),
4341  lower_band_width,upper_band_width,
4342  tol);
4343 
4344  }
4362  template <typename MatrixIterator1>
4363  inline
4364  bool is_tridiagonal(MatrixIterator1 A_up,
4365  MatrixIterator1 A_bot,
4366  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4367  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4368  {
4369  return slip::is_band_matrix(A_up,A_bot,1,1,tol);
4370  }
4371 
4372 
4388  template <typename Matrix>
4389  inline
4390  bool is_tridiagonal(const Matrix & A,
4391  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4392  = slip::epsilon<typename Matrix::value_type>())
4393  {
4394 
4395  return slip::is_tridiagonal(A.upper_left(),A.bottom_right(),tol);
4396 
4397  }
4398 
4416  template <typename MatrixIterator1>
4417  inline
4418  bool is_upper_bidiagonal(MatrixIterator1 A_up,
4419  MatrixIterator1 A_bot,
4420  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4421  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4422  {
4423  return slip::is_band_matrix(A_up,A_bot,0,1,tol);
4424  }
4440  template <typename Matrix>
4441  inline
4442  bool is_upper_bidiagonal(const Matrix & A,
4443  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4444  = slip::epsilon<typename Matrix::value_type>())
4445  {
4446 
4447  return slip::is_upper_bidiagonal(A.upper_left(),A.bottom_right(),tol);
4448 
4449  }
4450 
4468  template <typename MatrixIterator1>
4469  inline
4470  bool is_lower_bidiagonal(MatrixIterator1 A_up,
4471  MatrixIterator1 A_bot,
4472  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4473  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4474  {
4475  return slip::is_band_matrix(A_up,A_bot,1,0,tol);
4476  }
4477 
4493  template <typename Matrix>
4494  inline
4495  bool is_lower_bidiagonal(const Matrix & A,
4496  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4497  = slip::epsilon<typename Matrix::value_type>())
4498  {
4499 
4500  return slip::is_lower_bidiagonal(A.upper_left(),A.bottom_right(),tol);
4501 
4502  }
4519  template <typename MatrixIterator1>
4520  inline
4521  bool is_upper_hessenberg(MatrixIterator1 A_up,
4522  MatrixIterator1 A_bot,
4523  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4524  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4525  {
4526  return slip::is_band_matrix(A_up,A_bot,1,((A_bot-A_up)[0]-1),tol);
4527  }
4528 
4529 
4545  template <typename Matrix>
4546  inline
4547  bool is_upper_hessenberg(const Matrix & A,
4548  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4549  = slip::epsilon<typename Matrix::value_type>())
4550  {
4551 
4552  return slip::is_upper_hessenberg(A.upper_left(),A.bottom_right(),tol);
4553 
4554  }
4555 
4573  template <typename MatrixIterator1>
4574  inline
4575  bool is_lower_hessenberg(MatrixIterator1 A_up,
4576  MatrixIterator1 A_bot,
4577  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4578  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4579  {
4580  return slip::is_band_matrix(A_up,A_bot,((A_bot-A_up)[0]-1),1,tol);
4581  }
4582 
4583 
4599  template <typename Matrix>
4600  inline
4601  bool is_lower_hessenberg(const Matrix & A,
4602  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4603  = slip::epsilon<typename Matrix::value_type>())
4604  {
4605 
4606  return slip::is_lower_hessenberg(A.upper_left(),A.bottom_right(),tol);
4607 
4608  }
4609 
4626  template <typename MatrixIterator1>
4627  inline
4628  bool is_upper_triangular(MatrixIterator1 A_up,
4629  MatrixIterator1 A_bot,
4630  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4631  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4632  {
4633  return slip::is_band_matrix(A_up,A_bot,0,((A_bot-A_up)[0]-1),tol);
4634  }
4635 
4651  template <typename Matrix>
4652  inline
4653  bool is_upper_triangular(const Matrix & A,
4654  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4655  = slip::epsilon<typename Matrix::value_type>())
4656  {
4657 
4658  return slip::is_upper_triangular(A.upper_left(),A.bottom_right(),tol);
4659 
4660  }
4677  template <typename MatrixIterator1>
4678  inline
4679  bool is_lower_triangular(MatrixIterator1 A_up,
4680  MatrixIterator1 A_bot,
4681  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4682  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4683  {
4684  return slip::is_band_matrix(A_up,A_bot,((A_bot-A_up)[0]-1),0,tol);
4685  }
4701  template <typename Matrix>
4702  inline
4703  bool is_lower_triangular(const Matrix & A,
4704  const typename slip::lin_alg_traits<typename Matrix::value_type>::value_type tol
4705  = slip::epsilon<typename Matrix::value_type>())
4706  {
4707 
4708  return slip::is_lower_triangular(A.upper_left(),A.bottom_right(),tol);
4709 
4710  }
4728  template <typename MatrixIterator1>
4729  inline
4730  bool is_null_diagonal(MatrixIterator1 A_up,
4731  MatrixIterator1 A_bot,
4732  int diag_number,
4733  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4734  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4735 {
4736 
4737  int rows = int((A_bot - A_up)[0]);
4738  int cols = int((A_bot - A_up)[1]);
4739  int iterations = std::min(rows,cols) - std::abs(diag_number);
4740 
4741  slip::DPoint2d<int> d(1,1);
4742  int i = 0;
4743  int j = 0;
4744  if(diag_number >= 0)
4745  {
4746  i = 0;
4747  j = diag_number;
4748  }
4749  else
4750  {
4751  i = -diag_number;
4752  j = 0;
4753  }
4754  slip::DPoint2d<int> dstart(i,j);
4755  MatrixIterator1 A_up_start = A_up + dstart;
4756  for(int k = 0; k < iterations; ++k)
4757  {
4758  if(std::abs(*A_up_start) > tol)
4759  {
4760  return false;
4761  }
4762  A_up_start += d;
4763  }
4764  return true;
4765 }
4766 
4783  template <typename MatrixIterator1>
4784  inline
4785  slip::DENSE_MATRIX_TYPE analyse_matrix(MatrixIterator1 A_up,
4786  MatrixIterator1 A_bot,
4787  const typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator1>::value_type>::value_type tol
4788  = slip::epsilon<typename std::iterator_traits<MatrixIterator1>::value_type>())
4789  {
4790  int rows = int((A_bot - A_up)[0]);
4791  int cols = int((A_bot - A_up)[1]);
4792 
4793  int up_band_number = (cols - 1);
4794  int low_band_number = (rows - 1);
4795 
4796  int k = -low_band_number;
4797  while((k < 0) && slip::is_null_diagonal(A_up,A_bot,k,tol))
4798  {
4799  low_band_number--;
4800  ++k;
4801  }
4802  k = up_band_number;
4803  while((k > 0) && slip::is_null_diagonal(A_up,A_bot,k,tol))
4804  {
4805  up_band_number--;
4806  --k;
4807  }
4808 
4809  if((up_band_number == low_band_number) && (up_band_number = 0))
4810  {
4811  if(slip::is_null_diagonal(A_up,A_bot,0,tol))
4812  {
4813  return slip::DENSE_MATRIX_NULL;
4814  }
4815  else
4816  {
4817  return slip::DENSE_MATRIX_DIAGONAL;
4818  }
4819  }
4820  else if((up_band_number == (cols-1)) && (low_band_number == rows - 1))
4821  {
4822  return slip::DENSE_MATRIX_FULL;
4823  }
4824  else if((up_band_number == 1) && (low_band_number == 0))
4825  {
4826  return slip::DENSE_MATRIX_UPPER_BIDIAGONAL;
4827  }
4828  else if((up_band_number == 0) && (low_band_number == 1))
4829  {
4830  return slip::DENSE_MATRIX_LOWER_BIDIAGONAL;
4831  }
4832  else if((up_band_number == 1) && (low_band_number == 1))
4833  {
4834  return slip::DENSE_MATRIX_TRIDIAGONAL;
4835  }
4836  else if((up_band_number == (cols-1)) && (low_band_number == 0))
4837  {
4838  return slip::DENSE_MATRIX_UPPER_TRIANGULAR;
4839  }
4840  else if((up_band_number == 0) && (low_band_number == (rows-1)))
4841  {
4842  return slip::DENSE_MATRIX_LOWER_TRIANGULAR;
4843  }
4844  else if((up_band_number == (cols-1)) && (low_band_number == 1))
4845  {
4846  return slip::DENSE_MATRIX_UPPER_HESSENBERG;
4847  }
4848  else if((up_band_number == 1) && (low_band_number == (rows-1)))
4849  {
4850  return slip::DENSE_MATRIX_LOWER_HESSENBERG;
4851  }
4852  else
4853  {
4854  return slip::DENSE_MATRIX_BANDED;
4855  }
4856 
4857  }
4858 /* @} */
4859 
4860 
4862  /* @{ */
4863 
4864 
4875  template <typename MatrixIterator>
4876  inline
4877  int pivot_max(MatrixIterator M_up,MatrixIterator M_bot)
4878  {
4879  typedef typename std::iterator_traits<MatrixIterator>::value_type _T;
4880  typedef typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator>::value_type>::value_type norm_type;
4881  typename MatrixIterator::difference_type sizeM= M_bot - M_up;
4882 
4883  norm_type max = std::abs(_T());
4884  int ind = -1;
4885 
4886  while(M_up.x1() < M_bot.x1())
4887  {
4888  if(std::abs(*M_up) > max)
4889  {
4890  max = std::abs(*M_up);
4891  ind = M_up.x1();
4892  }
4893  M_up++;
4894  }
4895  return ind;
4896  }
4907 template <typename MatrixIterator>
4908 inline
4909 int partial_pivot(MatrixIterator M_up, MatrixIterator M_bot)
4910 {
4911  typedef typename std::iterator_traits<MatrixIterator>::value_type _T;
4912  typedef typename slip::lin_alg_traits<typename std::iterator_traits<MatrixIterator>::value_type>::value_type norm_type;
4913 
4914  norm_type max = std::abs(*M_up);
4915  int ind = 0;
4916  slip::DPoint2d<int> d10(1,0);
4917  MatrixIterator it = M_up + d10;
4918  int istop = int((M_bot - M_up)[0]);
4919  for(int i = 1; i < istop; ++i, it+=d10)
4920  {
4921  if(std::abs(*it) > max)
4922  {
4923  ind = i;
4924  max = std::abs(*it);
4925  }
4926  }
4927  if(max < std::abs(_T()))
4928  {
4929  ind = -1;
4930  }
4931  else
4932  {
4933  ind = ind + M_up.x1();
4934  }
4935  return ind;
4936 
4937 }
4938 /* @} */
4939 
4941  /* @{ */
4961  template <class Matrix1,class Matrix2>
4962  inline
4963  void inverse(const Matrix1 & M,
4964  const std::size_t nr1,
4965  const std::size_t nc1,
4966  Matrix2 & IM,
4967  const std::size_t nr2,
4968  const std::size_t nc2)
4969  {
4970  assert(nr1==nc1);
4971  assert(nr2==nc2);
4972  assert(nr1==nr2);
4973 
4974  size_t i,j,k;
4975 
4976  typedef typename Matrix2::value_type _T;
4977  _T tmpval;
4978 
4979 
4980  Array2d<typename Matrix2::value_type> Mcopy(nr1,nc1);
4981 
4982 
4983  // Copy M to the first half of the
4984  for (i=0;i<nr1;++i)
4985  {
4986  for (j=0;j<nc1;++j)
4987  {
4988  Mcopy[i][j]=M[i][j];
4989  }
4990  }
4991 
4992  // Identity matrix to the right
4993  std::fill(IM.begin(),IM.end(),_T(0));
4994  for (i=0;i<nr2;++i)
4995  {
4996  IM[i][i]=_T(1.0);
4997  }
4998 
4999 
5000  for (i=0;i<nr1;++i)
5001  {
5002  // M[i][i]==0 : find another row to add...
5003  if (Mcopy[i][i]==_T(0))
5004  {
5005  k=i+1;
5006  while ( k<nr1 && Mcopy[k][i]==_T(0) )
5007  {
5008  k++;
5009  }
5010  if ( k==nr1 )
5011  {
5012  // TODO : Need to raise exception
5013  std::cerr<< "Singular Matrix - can not invert it"<<std::endl;
5014  exit(1);
5015  }
5016  for ( j = i ; j < nr1;++j)
5017  {
5018  Mcopy[i][j]+= Mcopy[k][j] / Mcopy[k][i];
5019  }
5020  for ( j = 0 ; j < nr2;++j)
5021  {
5022  IM[i][j]+= IM[k][j] / Mcopy[k][i];
5023  }
5024 
5025  }
5026  // Scale the row to 1
5027  tmpval=Mcopy[i][i];
5028  for ( j = i ; j < nr1;++j)
5029  {
5030  Mcopy[i][j]/=tmpval;
5031  }
5032  for ( j = 0 ; j < nr1;++j)
5033  {
5034  IM[i][j]/=tmpval;
5035  }
5036  // Do elimination
5037  for ( k = 0 ; k < nr1;++k)
5038  {
5039  if (k != i)
5040  {
5041  tmpval=Mcopy[k][i];
5042  for ( j = i ; j < nr1;++j)
5043  {
5044  Mcopy[k][j]-=tmpval/Mcopy[i][i]*Mcopy[i][j];
5045  }
5046  for ( j = 0 ; j < nr2;++j)
5047  {
5048  IM[k][j]-=tmpval/Mcopy[i][i]*IM[i][j];
5049  }
5050  }
5051  }
5052  }
5053  }
5054 
5081  template <class Matrix1, class Matrix2>
5082  inline
5083  void inverse(const Matrix1 & M, Matrix2 & IM)
5084  {
5085  slip::inverse(M,M.rows(),M.cols(),IM,IM.rows(),IM.cols());
5086  }
5087 /* @} */
5088 
5089 
5091 /* @{ */
5162  template <class Matrix1, class Matrix2, class Vector>
5163  inline
5164  int lu(const Matrix1 & M,
5165  Matrix2 & LU,
5166  Vector & Indx)
5167  {
5168  assert(M.rows()==M.cols());
5169  assert(LU.rows()==LU.cols());
5170  assert(M.rows()==LU.rows());
5171  assert(M.rows()==Indx.size());
5172 
5173  std::size_t nr = M.rows();
5174  std::size_t nc = M.cols();
5175 
5176 
5177  typedef typename Matrix1::value_type value_type;
5178  typedef typename slip::lin_alg_traits<value_type>::value_type norm_type;
5179 
5180  slip::Array<norm_type> Vv(nr);
5181  int npivot = 1;
5182 
5183  for(std::size_t i = 0; i < nr; ++i)
5184  {
5185  for(std::size_t j = 0; j < nc; ++j)
5186  {
5187  LU[i][j] = M[i][j];
5188  }
5189  Indx[i] = i;
5190  }
5191  // Crout's algorithm without pivoting : working on LU
5192  for(std::size_t i = 0; i < nr; ++i)
5193  {
5194  norm_type max = std::abs(LU[i][0]);
5195  for(std::size_t j = 1; j < nc; ++j)
5196  {
5197  norm_type tmp2 =
5198  std::abs(LU[i][j]);
5199  if (tmp2 > max)
5200  {
5201  max = tmp2;
5202  }
5203  }
5204  if (max == typename slip::lin_alg_traits<value_type>::value_type(0))
5205  {
5206  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
5207  }
5208  Vv[i] = norm_type(1) / max;
5209  }
5210 
5211  for (std::size_t j = 0; j < nr; ++j)
5212  {
5213  for (std::size_t i = 0; i < j; ++i)
5214  {
5215  value_type sum = LU[i][j];
5216  for (std::size_t k = 0; k < i; ++k)
5217  {
5218  sum -= LU[i][k] * LU[k][j];
5219  }
5220  LU[i][j] = sum;
5221  }
5222  norm_type max = 0;
5223  std::size_t imax = 0;
5224  for (std::size_t i = j; i < nr; ++i)
5225  {
5226  value_type sum = LU[i][j];
5227  for (std::size_t k = 0; k < j; ++k)
5228  {
5229  sum -= LU[i][k] * LU[k][j];
5230  }
5231  LU[i][j] = sum;
5232  norm_type tmp = Vv[i] * std::abs(sum);
5233  if (tmp >= max)
5234  {
5235  max = tmp;
5236  imax = i;
5237  }
5238  }
5239  // Invert columns
5240  if (j != imax)
5241  {
5242  for (std::size_t k = 0; k < nr; ++k)
5243  {
5244  std::swap(LU[imax][k],LU[j][k]);
5245  }
5246  Vv[imax] = Vv[j];
5247 
5248  std::swap(Indx[j],Indx[imax]);
5249  npivot = - npivot;
5250  }
5251 
5252  //replace 0 values by epsilon
5253  if (LU[j][j] == value_type(0))
5254  {
5255  LU[j][j] = std::numeric_limits<value_type>::epsilon();
5256  }
5257  if (j != (nr-1))
5258  {
5259  value_type scal = value_type(1) / LU[j][j];
5260  for (std::size_t i =j+1; i < nr; ++i)
5261  {
5262  LU[i][j]*=scal;
5263  }
5264  }
5265  }
5266  return npivot;
5267 
5268  }
5269 
5270 
5351  template <class Matrix1, class Matrix2>
5352  inline
5353  int lu(const Matrix1 & M,
5354  Matrix2 & L,
5355  Matrix2 & U,
5356  Matrix2 & P)
5357  {
5358  assert(M.rows() == M.cols());
5359  assert(L.rows() == L.cols());
5360  assert(U.rows() == U.cols());
5361  assert(P.rows() == P.cols());
5362  assert(M.rows() == L.rows());
5363  assert(M.rows() == U.rows());
5364  assert(M.rows() == P.cols());
5365 
5366  std::size_t nr = M.rows();
5367  slip::Array<std::size_t> Indx(nr);
5368  //decompostion LU
5369  int sign = slip::lu(M,U,Indx);
5370  // Generate L and P
5371  for (std::size_t i = 0; i < nr; ++i)
5372  {
5373  for (std::size_t j=0; j < i; ++j)
5374  {
5375  L[i][j] = U[i][j];
5376  U[i][j] = 0;
5377  P[i][j] = 0;
5378  }
5379  }
5380  for (std::size_t i = 0; i < nr; ++i)
5381  {
5382  L[i][i] = 1;
5383  P[Indx[i]][i] = 1;
5384  }
5385  return sign;
5386  }
5387 
5388 
5419  template <class Matrix, class Vector1, class Vector2, class Vector3>
5420  inline
5421  void lu_solve(const Matrix LU,
5422  const std::size_t nr,
5423  const std::size_t nc,
5424  const Vector1 & Indx,
5425  const std::size_t nv1,
5426  const Vector2 & B,
5427  const std::size_t nv2,
5428  Vector3 & X,
5429  const std::size_t nv3)
5430  {
5431  assert(nr==nc);
5432  assert(nv1==nv2);
5433  assert(nv2==nv3);
5434  assert(nr==nv1);
5435 
5436  // Solving LY=B
5437  for (std::size_t i = 0; i < nr; ++i)
5438  {
5439  // Permut value
5440  X[i] = B[Indx[i]];
5441 
5442  for (std::size_t j = 0; j < i; ++j)
5443  {
5444  X[i] -= LU[i][j] * X[j];
5445 
5446  }
5447  }
5448  // Solving UX=Y : std::size_t always positive so check if <nr for loop ending
5449  for (std::size_t i = nr-1; i < nr; i--)
5450  {
5451  for (std::size_t j=i+1;j<nr;++j)
5452  {
5453  X[i] -= LU[i][j] * X[j];
5454  }
5455  X[i] /= LU[i][i];
5456  }
5457 
5458  }
5459 
5493  template <class Matrix, class Vector1, class Vector2>
5494  inline
5495  void lu_solve(const Matrix& A,
5496  Vector1& X,
5497  const Vector2& B)
5498  {
5499  std::size_t nr = A.rows();
5500  std::size_t nc = A.cols();
5501 
5502  Matrix LU(nr,nc);
5503  slip::Array<std::size_t> Indx(nr);
5504  slip::lu(A,LU,Indx);
5505  // Solving LY=B
5506  for (std::size_t i = 0; i < nr; ++i)
5507  {
5508  // Permut value
5509  X[i] = B[Indx[i]];
5510 
5511  for (std::size_t j = 0; j < i; ++j)
5512  {
5513  X[i] -= LU[i][j] * X[j];
5514  }
5515  }
5516  // Solving UX=Y:
5517  for (std::size_t i = nr-1; i < nr; i--)
5518  {
5519  for (std::size_t j = i+1;j < nr; ++j)
5520  {
5521  X[i] -= LU[i][j] * X[j];
5522  }
5523  X[i] /= LU[i][i];
5524  }
5525  }
5526 
5527 
5556  template <class Matrix1,class Matrix2>
5557  inline
5558  void lu_inv(const Matrix1 & M,
5559  Matrix2 & IM)
5560  {
5561  assert(M.rows() == M.cols());
5562  assert(IM.rows() == IM.cols());
5563  assert(M.rows() == IM.rows());
5564 
5565  std::size_t nr1 = M.rows();
5566  std::size_t nc1 = M.cols();
5567 
5568  typedef typename Matrix1::value_type value_type;
5569 
5570  slip::Array<std::size_t> permut(nr1);
5571  slip::Array2d<value_type> LU(nr1,nc1);
5572  slip::Array<value_type> B(nr1);
5573  slip::Array<value_type> X(nr1);
5574 
5575 
5576 
5577  lu(M,LU,permut);
5578 
5579  for (std::size_t j = 0; j < nc1; ++j)
5580  {
5581  for (std::size_t i = 0; i < nr1; ++i)
5582  {
5583  B[i] = value_type(0.0);
5584  }
5585  B[j] = value_type(1.0);
5586  lu_solve(LU,nr1,nc1,permut,nr1,B,nr1,X,nr1);
5587  for (std::size_t i = 0;i < nr1; ++i)
5588  {
5589  IM[i][j] = X[i];
5590  }
5591  }
5592  }
5593 
5594 
5595 
5616  template <class Matrix1>
5617  inline
5618  typename Matrix1::value_type
5619  lu_det(const Matrix1& M)
5620  {
5621  assert(M.rows() == M.cols());
5622  std::size_t nr1 = M.rows();
5623  std::size_t nc1 = M.cols();
5624 
5625  typedef typename Matrix1::value_type value_type;
5626  value_type det = 1.0;
5627  int sign = 0;
5628  slip::Array2d<value_type> LU(nr1,nc1);
5629  slip::Array<std::size_t> permut(nr1);
5630  sign = lu(M,LU,permut);
5631  for (std::size_t i=0;i<nr1;++i)
5632  {
5633  det *= LU[i][i];
5634  }
5635 
5636  if(sign < 0)
5637  {
5638  det = -det;
5639  }
5640  return det;
5641  }
5642 
5643 
5644 /* @} */
5645 
5647  /* @{ */
5648 
5720  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
5721  inline
5722  void upper_triangular_solve(MatrixIterator U_up,
5723  MatrixIterator U_bot,
5724  RandomAccessIterator1 X_first,
5725  RandomAccessIterator1 X_last,
5726  RandomAccessIterator2 B_first,
5727  RandomAccessIterator2 B_last,
5728  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision)
5729  {
5730  assert((U_bot - U_up)[0] == (U_bot - U_up)[1]);
5731  assert((U_bot - U_up)[1] == (X_last - X_first));
5732  assert((X_last - X_first) == (B_last - B_first));
5733 
5734 
5735  typedef typename MatrixIterator::value_type _T;
5736 
5737  int U_rows_m1 = int((U_bot - U_up)[0]) - 1;
5738  int U_rows_m2 = U_rows_m1 - 1;
5739  slip::DPoint2d<int> d(U_rows_m1,U_rows_m1);
5740 
5741  //bottom right value inside the box
5742  _T botr = U_up[d];
5743  if(std::abs(botr) < precision)
5744  {
5745 
5746  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
5747  }
5748  //computes the last value of X
5749  *(--X_last) = *(--B_last) / botr;
5750  --X_last;
5751  --B_last;
5752  for(int i = U_rows_m2; i >= 0 ; i--, X_last--, B_last--)
5753  {
5754  d.set_coord(i,i);
5755  _T uii = U_up[d];
5756  _T inv_uii = _T(1.0)/uii;
5757  if(std::abs(uii) < precision)
5758  {
5759  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
5760  }
5761  int istart = (i+1);
5762  _T sum = std::inner_product(U_up.row_begin(i) + istart,U_up.row_end(i),
5763  X_first + istart,_T());
5764  *X_last = (*B_last - sum) * inv_uii;
5765  }
5766 
5767 
5768  }
5769 
5770 
5834  template <class Matrix,class Vector1,class Vector2>
5835  inline
5836  void upper_triangular_solve(const Matrix & U,
5837  Vector1 & X,
5838  const Vector2 & B,
5839  typename slip::lin_alg_traits<typename Matrix::value_type>::value_type precision)
5840  {
5841  slip::upper_triangular_solve(U.upper_left(),U.bottom_right(),
5842  X.begin(),X.end(),
5843  B.begin(),B.end(),
5844  precision);
5845  }
5846 
5892  template <typename MatrixIterator1,
5893  typename MatrixIterator2>
5894  inline
5895  void upper_triangular_inv(MatrixIterator1 A_up,
5896  MatrixIterator1 A_bot,
5897  MatrixIterator2 Ainv_up,
5898  MatrixIterator2 Ainv_bot,
5899  typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type precision = typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type(1.0E-6))
5900  {
5901  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
5902  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
5903  assert((A_bot - A_up)[0] == (Ainv_bot - Ainv_up)[0]);
5904  assert((A_bot - A_up)[1] == (Ainv_bot - Ainv_up)[1]);
5905 
5906 
5907  typedef typename MatrixIterator1::value_type value_type;
5908  std::size_t n = std::size_t((A_bot - A_up)[0]);
5909  //init Ainv with the identity matrix
5910  slip::identity(Ainv_up,Ainv_bot);
5911 
5912 
5913  //solve the n systems
5914  for(std::size_t k = 0; k < n; ++k)
5915  {
5916  //resolution of LY = B
5917  slip::upper_triangular_solve(A_up,
5918  A_bot,
5919  Ainv_up.col_begin(k),
5920  Ainv_up.col_end(k),
5921  Ainv_up.col_begin(k),
5922  Ainv_up.col_end(k),
5923  precision);
5924 
5925  }
5926 
5927 
5928  }
5967  template <typename Matrix1,
5968  typename Matrix2>
5969  inline
5970  void upper_triangular_inv(Matrix1 A,
5971  Matrix2 Ainv,
5972  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type(1.0E-6))
5973  {
5974  slip::upper_triangular_inv(A.upper_left(),A.bottom_right(),
5975  Ainv.upper_left(),Ainv.bottom_right(),
5976  precision);
5977  }
6046  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
6047  inline
6048  void lower_triangular_solve(MatrixIterator L_up,
6049  MatrixIterator L_bot,
6050  RandomAccessIterator1 X_first,
6051  RandomAccessIterator1 X_last,
6052  RandomAccessIterator2 B_first,
6053  RandomAccessIterator2 B_last,
6054  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision)
6055  {
6056  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
6057  assert((L_bot - L_up)[1] == (X_last - X_first));
6058  assert((X_last - X_first) == (B_last - B_first));
6059 
6060 
6061  typedef typename MatrixIterator::value_type _T;
6062 
6063  RandomAccessIterator1 X_first_tmp = X_first;
6064  int rows = int((L_bot - L_up)[0]);
6065  _T L00 = *L_up;
6066  if(std::abs(L00) < precision)
6067  {
6068  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6069  }
6070  slip::DPoint2d<int> d(1,1);
6071  *(X_first++) = *(B_first++) / L00;
6072  for(int i = 1; i < rows ; ++i, ++X_first, ++B_first)
6073  {
6074  d.set_coord(i,i);
6075  _T Lii = L_up[d];
6076  _T inv_Lii = _T(1.0)/Lii;
6077  if(std::abs(Lii) < precision)
6078  {
6079  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6080  }
6081  _T sum = std::inner_product(L_up.row_begin(i),L_up.row_begin(i)+i,X_first_tmp,_T());
6082 
6083  *X_first = (*B_first - sum) * inv_Lii;
6084  }
6085 
6086 
6087  }
6088 
6151  template <class Matrix,class Vector1,class Vector2>
6152  inline
6153  void lower_triangular_solve(const Matrix & L,
6154  Vector1 & X,
6155  const Vector2 & B,
6156  typename slip::lin_alg_traits<typename Matrix::value_type>::value_type precision)
6157  {
6158  slip::lower_triangular_solve(L.upper_left(),L.bottom_right(),
6159  X.begin(),X.end(),
6160  B.begin(),B.end(),
6161  precision);
6162  }
6163 
6164 
6209  template <typename MatrixIterator1,
6210  typename MatrixIterator2>
6211  inline
6212  void lower_triangular_inv(MatrixIterator1 A_up,
6213  MatrixIterator1 A_bot,
6214  MatrixIterator2 Ainv_up,
6215  MatrixIterator2 Ainv_bot,
6216  typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type precision = typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type(1.0E-6))
6217  {
6218  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
6219  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
6220  assert((A_bot - A_up)[0] == (Ainv_bot - Ainv_up)[0]);
6221  assert((A_bot - A_up)[1] == (Ainv_bot - Ainv_up)[1]);
6222 
6223 
6224  typedef typename MatrixIterator1::value_type value_type;
6225  std::size_t n = std::size_t((A_bot - A_up)[0]);
6226  //init Ainv with the identity matrix
6227  slip::identity(Ainv_up,Ainv_bot);
6228 
6229 
6230  //solve the n systems
6231  for(std::size_t k = 0; k < n; ++k)
6232  {
6233  //resolution of LY = B
6234  slip::lower_triangular_solve(A_up,
6235  A_bot,
6236  Ainv_up.col_begin(k),
6237  Ainv_up.col_end(k),
6238  Ainv_up.col_begin(k),
6239  Ainv_up.col_end(k),
6240  precision);
6241 
6242  }
6243  }
6244 
6282  template <typename Matrix1,
6283  typename Matrix2>
6284  inline
6285  void lower_triangular_inv(const Matrix1& A,
6286  Matrix2& Ainv,
6287  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type(1.0E-6))
6288  {
6289  slip::lower_triangular_inv(A.upper_left(),A.bottom_right(),
6290  Ainv.upper_left(),Ainv.bottom_right(),
6291  precision);
6292  }
6361  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
6362  inline
6363  void unit_upper_triangular_solve(MatrixIterator U_up,
6364  MatrixIterator U_bot,
6365  RandomAccessIterator1 X_first,
6366  RandomAccessIterator1 X_last,
6367  RandomAccessIterator2 B_first,
6368  RandomAccessIterator2 B_last)
6369  {
6370  assert((U_bot - U_up)[0] == (U_bot - U_up)[1]);
6371  assert((U_bot - U_up)[1] == (X_last - X_first));
6372  assert((X_last - X_first) == (B_last - B_first));
6373 
6374 
6375  typedef typename MatrixIterator::value_type _T;
6376 
6377  int U_rows_m1 = int((U_bot - U_up)[0]) - 1;
6378  int U_rows_m2 = U_rows_m1 - 1;
6379  slip::DPoint2d<int> d(U_rows_m1,U_rows_m1);
6380 
6381 
6382  //computes the last value of X
6383  *(--X_last) = *(--B_last);
6384  --X_last;
6385  --B_last;
6386  for(int i = U_rows_m2; i >= 0 ; i--, X_last--, B_last--)
6387  {
6388  d.set_coord(i,i);
6389 
6390  int istart = (i+1);
6391  _T sum = std::inner_product(U_up.row_begin(i) + istart,
6392  U_up.row_end(i),
6393  X_first + istart,_T());
6394  *X_last = (*B_last - sum) ;
6395  }
6396  }
6397 
6459  template <class Matrix,class Vector1,class Vector2>
6460  inline
6461  void unit_upper_triangular_solve(const Matrix & U,
6462  Vector1 & X,
6463  const Vector2 & B)
6464  {
6465 
6466  slip::unit_upper_triangular_solve(U.upper_left(),U.bottom_right(),
6467  X.begin(),X.end(),
6468  B.begin(),B.end());
6469  }
6470 
6471 
6537  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
6538  inline
6539  void unit_lower_triangular_solve(MatrixIterator L_up,
6540  MatrixIterator L_bot,
6541  RandomAccessIterator1 X_first,
6542  RandomAccessIterator1 X_last,
6543  RandomAccessIterator2 B_first,
6544  RandomAccessIterator2 B_last)
6545  {
6546  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
6547  assert((L_bot - L_up)[1] == (X_last - X_first));
6548  assert((X_last - X_first) == (B_last - B_first));
6549 
6550 
6551  typedef typename MatrixIterator::value_type _T;
6552 
6553  RandomAccessIterator1 X_first_tmp = X_first;
6554  int rows = int((L_bot - L_up)[0]);
6555 
6556  slip::DPoint2d<int> d(1,1);
6557  *(X_first++) = *(B_first++);
6558  for(int i = 1; i < rows ; ++i, ++X_first, ++B_first)
6559  {
6560  d.set_coord(i,i);
6561 
6562  _T sum = std::inner_product(L_up.row_begin(i),L_up.row_begin(i)+i,X_first_tmp,_T());
6563 
6564  *X_first = (*B_first - sum);
6565  }
6566  }
6567 
6568 
6629  template <class Matrix,class Vector1,class Vector2>
6630  inline
6631  void unit_lower_triangular_solve(const Matrix & L,
6632  Vector1 & X,
6633  const Vector2 & B)
6634  {
6635  slip::unit_lower_triangular_solve(L.upper_left(),L.bottom_right(),
6636  X.begin(),X.end(),
6637  B.begin(),B.end());
6638  }
6639 
6640 /* @} */
6641 
6643  /* @{ */
6712  template <typename RandomAccessIterator1,
6713  typename RandomAccessIterator2,
6714  typename RandomAccessIterator3>
6715  inline
6716  void diagonal_solve(RandomAccessIterator1 first,
6717  RandomAccessIterator1 last,
6718  RandomAccessIterator2 X_first,
6719  RandomAccessIterator2 X_last,
6720  RandomAccessIterator3 B_first,
6721  RandomAccessIterator3 B_last,
6722  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision)
6723  {
6724 
6725  assert((last - first) == (X_last - X_first));
6726  assert((B_last - B_first) == (X_last - X_first));
6727 
6728 
6729  //Computes X
6730  for(; first != last; ++first, ++X_first, ++B_first)
6731  {
6732  if(std::abs(*first) < precision)
6733  {
6734  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6735 
6736  }
6737  *X_first = *B_first / *first;
6738  }
6739  }
6740 /* @} */
6741 
6743  /* @{ */
6744 
6770 template <typename MatrixIterator1,
6771  typename RandomAccessIterator>
6772  inline
6773  int band_lu(MatrixIterator1 A_up,
6774  MatrixIterator1 A_bot,
6775  const int p,
6776  const int q,
6777  RandomAccessIterator Ind_first,
6778  RandomAccessIterator Ind_last,
6779  typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type precision )
6780  {
6781  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
6782  assert((A_bot - A_up)[1] == (Ind_last - Ind_first));
6783  assert(p < int((A_bot - A_up)[0]));
6784  assert(q < int((A_bot - A_up)[0]));
6785 
6786  typedef typename MatrixIterator1::value_type value_type;
6787  typedef typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type real_type;
6788  int A_rows = int((A_bot - A_up)[0]);
6789  int A_rows_m1 = A_rows - 1;
6790  slip::DPoint2d<int> d(0,0);
6791  slip::DPoint2d<int> d2(0,0);
6792  slip::DPoint2d<int> dstart(0,0);
6793  slip::DPoint2d<int> dstop(0,0);
6794  int nb_pivot = 1;
6795  slip::iota(Ind_first,Ind_last,0);
6796 
6797  //value_type proj = value_type();
6798  for(int k = 0; k < A_rows_m1; ++k)
6799  {
6800  d.set_coord(k,k);
6801  int istart = k + 1;
6802  int istop = std::min(k+p,A_rows_m1) + 1;
6803  int jstart = istart;
6804  int jstop = std::min(k+p+q,A_rows_m1) + 1;
6805  dstart.set_coord(k,k);
6806  dstop.set_coord(istop,jstop);
6807 
6808  int pivot = slip::partial_pivot(A_up + dstart,A_up + dstop);
6809  if(pivot == -1)
6810  {
6811  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6812 
6813  }
6814 
6815  if((pivot > 0) && (pivot != k))
6816  {
6817  std::swap_ranges(A_up.row_begin(k),A_up.row_end(k),
6818  A_up.row_begin(pivot));
6819  nb_pivot -= 1;
6820  std::swap(Ind_first[k],Ind_first[pivot]);
6821  }
6822  value_type Akk = A_up[d];
6823  if(std::abs(Akk) < precision)
6824  {
6825  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6826  }
6827  value_type inv_Akk = value_type(1) / Akk;
6828  slip::vector_scalar_multiplies(A_up.col_begin(k)+istart,
6829  A_up.col_begin(k)+istop,
6830  inv_Akk,
6831  A_up.col_begin(k)+istart);
6832 
6833  for(int i = istart; i < istop; ++i)
6834  {
6835  d2.set_coord(i,k);
6836  slip::reduction(A_up.row_begin(i)+jstart,
6837  A_up.row_begin(i)+jstop,
6838  A_up.row_begin(k)+jstart,
6839  -A_up[d2]);
6840  }
6841  }
6842 
6843  d.set_coord(A_rows_m1,A_rows_m1);
6844  if(std::abs(A_up[d]) < precision)
6845  {
6846  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
6847 
6848  }
6849  return nb_pivot;
6850 
6851  }
6852 
6900 template <class Matrix1, class Matrix2, class Matrix3, class Matrix4>
6901  inline
6902  int band_lu(const Matrix1 & M,
6903  const int p,
6904  const int q,
6905  Matrix2 & L,
6906  Matrix3 & U,
6907  Matrix4 & P,
6908  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type precision )
6909  {
6910  assert(M.rows() == M.cols());
6911  assert(L.rows() == L.cols());
6912  assert(U.rows() == U.cols());
6913  assert(P.rows() == P.cols());
6914  assert(M.rows() == L.rows());
6915  assert(M.rows() == U.rows());
6916  assert(M.rows() == P.cols());
6917 
6918  std::size_t nr = M.rows();
6919  slip::Array<std::size_t> Indx(nr);
6920 
6921  std::copy(M.begin(),M.end(),U.begin());
6922  //decompostion LU
6923  int sign = slip::band_lu(U.upper_left(),
6924  U.bottom_right(),
6925  p,q,
6926  Indx.begin(),
6927  Indx.end(),
6928  precision);
6929  // Generate L and P
6930  for (std::size_t i = 0; i < nr; ++i)
6931  {
6932  for (std::size_t j=0; j < i; ++j)
6933  {
6934  L[i][j] = U[i][j];
6935  U[i][j] = 0;
6936  P[i][j] = 0;
6937  }
6938  }
6939  for (std::size_t i = 0; i < nr; ++i)
6940  {
6941  L[i][i] = 1;
6942  P[Indx[i]][i] = 1;
6943  }
6944  return sign;
6945  }
7023 template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
7024  inline
7025  void upper_band_solve(MatrixIterator U_up,
7026  MatrixIterator U_bot,
7027  int q,
7028  RandomAccessIterator1 X_first,
7029  RandomAccessIterator1 X_last,
7030  RandomAccessIterator2 B_first,
7031  RandomAccessIterator2 B_last,
7032  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision)
7033  {
7034  assert((U_bot - U_up)[0] == (U_bot - U_up)[1]);
7035  assert((U_bot - U_up)[1] == (X_last - X_first));
7036  assert((X_last - X_first) == (B_last - B_first));
7037  assert( q < int((U_bot - U_up)[0]));
7038 
7039 
7040  typedef typename MatrixIterator::value_type _T;
7041 
7042  int U_rows_m1 = int((U_bot - U_up)[0]) - 1;
7043  int U_rows_m2 = U_rows_m1 - 1;
7044  slip::DPoint2d<int> d(U_rows_m1,U_rows_m1);
7045 
7046  //bottom right value inside the box
7047  _T botr = U_up[d];
7048  if(std::abs(botr) < precision)
7049  {
7050 
7051  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7052  }
7053  //computes the last value of X
7054  *(--X_last) = *(--B_last) / botr;
7055  --X_last;
7056  --B_last;
7057  for(int i = U_rows_m2; i >= 0 ; i--, X_last--, B_last--)
7058  {
7059  d.set_coord(i,i);
7060  _T uii = U_up[d];
7061  _T inv_uii = _T(1.0)/uii;
7062  if(std::abs(uii) < precision)
7063  {
7064  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7065  }
7066  int istart = (i+1);
7067  int istop = std::min(i+q,U_rows_m1) + 1;
7068  _T sum = std::inner_product(U_up.row_begin(i) + istart,U_up.row_begin(i) + istop,X_first + istart,_T());
7069  *X_last = (*B_last - sum) * inv_uii;
7070  }
7071 
7072 
7073  }
7152  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
7153  inline
7154  void unit_upper_band_solve(MatrixIterator U_up,
7155  MatrixIterator U_bot,
7156  int q,
7157  RandomAccessIterator1 X_first,
7158  RandomAccessIterator1 X_last,
7159  RandomAccessIterator2 B_first,
7160  RandomAccessIterator2 B_last,
7161  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision)
7162  {
7163  assert((U_bot - U_up)[0] == (U_bot - U_up)[1]);
7164  assert((U_bot - U_up)[1] == (X_last - X_first));
7165  assert((X_last - X_first) == (B_last - B_first));
7166  assert( q < int((U_bot - U_up)[0]));
7167 
7168 
7169  typedef typename MatrixIterator::value_type _T;
7170 
7171  int U_rows_m1 = int((U_bot - U_up)[0]) - 1;
7172  int U_rows_m2 = U_rows_m1 - 1;
7173  slip::DPoint2d<int> d(U_rows_m1,U_rows_m1);
7174 
7175 
7176  //computes the last value of X
7177  *(--X_last) = *(--B_last);
7178  --X_last;
7179  --B_last;
7180  for(int i = U_rows_m2; i >= 0 ; i--, X_last--, B_last--)
7181  {
7182  d.set_coord(i,i);
7183  int istart = (i+1);
7184  int istop = std::min(i+q,U_rows_m1) + 1;
7185  _T sum = std::inner_product(U_up.row_begin(i) + istart,U_up.row_begin(i) + istop,X_first + istart,_T());
7186  *X_last = (*B_last - sum);
7187  }
7188 
7189 
7190  }
7191 
7192 
7271  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
7272  inline
7273  void lower_band_solve(MatrixIterator L_up,
7274  MatrixIterator L_bot,
7275  int p,
7276  RandomAccessIterator1 X_first,
7277  RandomAccessIterator1 X_last,
7278  RandomAccessIterator2 B_first,
7279  RandomAccessIterator2 B_last,
7280  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision)
7281  {
7282  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
7283  assert((L_bot - L_up)[1] == (X_last - X_first));
7284  assert((X_last - X_first) == (B_last - B_first));
7285  assert( p < int((L_bot - L_up)[0]));
7286 
7287 
7288  typedef typename MatrixIterator::value_type _T;
7289 
7290  RandomAccessIterator1 X_first_tmp = X_first;
7291  int rows = int((L_bot - L_up)[0]);
7292  _T L00 = *L_up;
7293  if(std::abs(L00) < precision)
7294  {
7295  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7296  }
7297  slip::DPoint2d<int> d(1,1);
7298  *(X_first++) = *(B_first++) / L00;
7299  for(int i = 1; i < rows ; ++i, ++X_first, ++B_first)
7300  {
7301  d.set_coord(i,i);
7302  _T Lii = L_up[d];
7303  _T inv_Lii = _T(1.0)/Lii;
7304  int istart = std::max(0,i-p);
7305  if(std::abs(Lii) < precision)
7306  {
7307  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7308  }
7309 
7310  _T sum = std::inner_product(L_up.row_begin(i)+istart,L_up.row_begin(i)+i,X_first_tmp+istart,_T());
7311 
7312  *X_first = (*B_first - sum) * inv_Lii;
7313  }
7314 
7315 
7316 
7317  }
7318 
7319 
7398  template <typename MatrixIterator,typename RandomAccessIterator1,typename RandomAccessIterator2>
7399  inline
7400  void unit_lower_band_solve(MatrixIterator L_up,
7401  MatrixIterator L_bot,
7402  int p,
7403  RandomAccessIterator1 X_first,
7404  RandomAccessIterator1 X_last,
7405  RandomAccessIterator2 B_first,
7406  RandomAccessIterator2 B_last)
7407  {
7408  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
7409  assert((L_bot - L_up)[1] == (X_last - X_first));
7410  assert((X_last - X_first) == (B_last - B_first));
7411  assert( p < int((L_bot - L_up)[0]));
7412 
7413 
7414  typedef typename MatrixIterator::value_type _T;
7415 
7416  RandomAccessIterator1 X_first_tmp = X_first;
7417  int rows = int((L_bot - L_up)[0]);
7418 
7419  slip::DPoint2d<int> d(1,1);
7420  *(X_first++) = *(B_first++);
7421  for(int i = 1; i < rows ; ++i, ++X_first, ++B_first)
7422  {
7423  d.set_coord(i,i);
7424 
7425  int istart = std::max(0,i-p);
7426 
7427  _T sum = std::inner_product(L_up.row_begin(i)+istart,L_up.row_begin(i)+i,X_first_tmp+istart,_T());
7428 
7429  *X_first = (*B_first - sum);
7430  }
7431 
7432 
7433  }
7434 
7435 
7436 //TO FIX!!!
7437 /*
7438  template <typename MatrixIterator,
7439  typename RandomAccessIterator1,
7440  typename RandomAccessIterator2>
7441  inline
7442  void band_lu_solve(MatrixIterator A_up,
7443  MatrixIterator A_bot,
7444  const int p,
7445  const int q,
7446  RandomAccessIterator1 X_first,
7447  RandomAccessIterator1 X_last,
7448  RandomAccessIterator2 B_first,
7449  RandomAccessIterator2 B_last,
7450  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision = typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type(1.0E-6))
7451  {
7452  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
7453  assert((A_bot - A_up)[1] == (X_last - X_first));
7454  assert((X_last - X_first) == (B_last - B_first));
7455 
7456  typedef typename std::iterator_traits<MatrixIterator>::value_type value_type;
7457  typedef typename MatrixIterator::size_type size_type;
7458  size_type n = size_type((A_bot - A_up)[0]);
7459 
7460  //decomposition
7461  slip::Array2d<value_type> LU(n,n);
7462  slip::Array<int> Indx(n);
7463  std::copy(A_up,A_bot,LU.begin());
7464  int sign = slip::band_lu(LU.upper_left(),
7465  LU.bottom_right(),
7466  p,q,
7467  Indx.begin(),
7468  Indx.end(),
7469  precision);
7470  std::cout<<"Indx = "<<Indx<<std::endl;
7471  slip::Array<typename std::iterator_traits<RandomAccessIterator2>::value_type > B2(n);
7472  //swap B according to Indx
7473  for(int i = 0; i < n; ++i)
7474  {
7475  //std::swap(B_first[i],B_first[Indx[i]]);
7476  //value_type tmp = B_first[i];
7477  B2[i] = B_first[Indx[i]];
7478  //B_first[Indx[i]] = tmp;
7479  }
7480  // std::copy(B_first,B_last,std::ostream_iterator<value_type>(std::cout," "));
7481 // std::cout<<std::endl;
7482  std::cout<<"B2 = "<<B2<<std::endl;
7483  slip::Array<value_type> Y(n);
7484  //resolution of LY = B
7485  // slip::unit_lower_band_solve(LU.upper_left(),LU.bottom_right(),
7486 // p+q,
7487 // Y.begin(),Y.end(),
7488 // B_first,B_last,precision);
7489  // slip::lower_triangular_solve(LU.upper_left(),LU.bottom_right(),
7490 // Y.begin(),Y.end(),
7491 // B_first,B_last,precision);
7492  slip::lower_triangular_solve(LU.upper_left(),LU.bottom_right(),
7493  Y.begin(),Y.end(),
7494  B2.begin(),B2.end(),precision);
7495  std::cout<<"Y = "<<Y<<std::endl;
7496  //resolution of UX = Y
7497  // slip::upper_band_solve(LU.upper_left(),LU.bottom_right(),
7498 // p+q,
7499 // X_first,X_last,
7500 // Y.begin(),Y.end(),precision);
7501  slip::upper_triangular_solve(LU.upper_left(),LU.bottom_right(),
7502  X_first,X_last,
7503  Y.begin(),Y.end(),precision);
7504 
7505  }
7506 
7507 */
7508 
7509 /* @} */
7510 
7512  /* @{ */
7513 
7564  template <class Matrix,class Vector1,class Vector2>
7565  inline
7566  bool gauss_solve(const Matrix & M,
7567  Vector1 & X,
7568  const Vector2 & B,
7569  typename slip::lin_alg_traits<typename Matrix::value_type>::value_type precision)
7570  {
7571 
7572  typedef typename Matrix::value_type _T;
7573  int M_rows = int(M.rows());
7574  int M_cols = int(M.cols());
7575  int M_rows_m1 = M_rows - 1;
7576  int M_cols_m1 = M_cols - 1;
7577  //int M_rows_p1 = M_rows + 1;
7578  int M_cols_p1 = M_cols + 1;
7579 
7580  //We first build S = [M|B]
7581  Matrix S(M_rows, M_cols_p1);
7582  slip::Box2d<int> bm(0,0,M_rows_m1,M_cols_m1);
7583  std::copy(M.upper_left(),M.bottom_right(),S.upper_left(bm));
7584  std::copy(B.begin(),B.end(),S.col_begin(M.cols()));
7585  int nbr = 0;
7586 
7587  //triangularization of S
7588  while(nbr < M_rows)
7589  {
7590  bm.set_coord(nbr,nbr,M_rows_m1,M_cols);
7591  int nbpivot = slip::pivot_max(S.upper_left(bm),S.bottom_right(bm));
7592  if(nbpivot == -1)
7593  {
7594  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7595 
7596  }
7597  if(std::abs(S[nbpivot][nbr]) < precision)
7598  {
7599  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7600 
7601  }
7602  if(nbpivot > 0)
7603  {
7604  std::swap_ranges(S.row_begin(nbr),S.row_end(nbr),S.row_begin(nbpivot));
7605  }
7606  for(int i = nbr+1; i < M_rows; ++i)
7607  {
7608  _T u = - S[i][nbr] / S[nbr][nbr];
7609  slip::reduction(S.row_begin(i),S.row_end(i),S.row_begin(nbr),u);
7610  }
7611  nbr++;
7612  }
7613 
7614  //Calculation of X
7615  if(std::abs(S[M_rows_m1][M_cols_m1]) < precision)
7616  {
7617  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7618 
7619  }
7620 
7621  bm.set_coord(0,0,M_rows_m1,M_cols_m1);
7622  slip::upper_triangular_solve(S.upper_left(bm),S.bottom_right(bm),
7623  X.begin(),X.end(),
7624  S.col_begin(M_cols),S.col_end(M_cols),
7625  precision);
7626  return true;
7627 
7628  }
7629 
7630 
7631 
7683  template <class Matrix,
7684  class Vector1,
7685  class Vector2>
7686  inline
7687  bool gauss_solve_with_filling(const Matrix & M,
7688  Vector1 & X,
7689  const Vector2 & B,
7690  typename slip::lin_alg_traits<typename Matrix::value_type>::value_type precision)
7691  {
7692  typedef typename Matrix::value_type _T;
7693  int M_rows = int(M.rows());
7694  int M_cols = int(M.cols());
7695  int M_rows_m1 = M_rows - 1;
7696  int M_cols_m1 = M_cols - 1;
7697  //int M_rows_p1 = M_rows + 1;
7698  int M_cols_p1 = M_cols + 1;
7699 
7700  //We first build S = [M|B]
7701  Matrix S(M_rows, M_cols_p1);
7702  slip::Box2d<int> bm(0,0,M_rows_m1,M_cols_m1);
7703  std::copy(M.upper_left(),M.bottom_right(),S.upper_left(bm));
7704  std::copy(B.begin(),B.end(),S.col_begin(M.cols()));
7705  int nbr = 0;
7706  bool ret = true;
7707 
7708  //triangularization of S
7709  while(nbr < M_rows)
7710  {
7711  bm.set_coord(nbr,nbr,M_rows_m1,M_cols);
7712  int nbpivot = slip::pivot_max(S.upper_left(bm),S.bottom_right(bm));
7713  if(nbpivot == -1)
7714  {
7715  ret = false;
7716  *(S.upper_left(bm)) = precision;
7717  nbpivot = S.upper_left(bm).x1();
7718 
7719  }
7720  if(std::abs(S[nbpivot][nbr]) < precision)
7721  {
7722  ret = false;
7723  S[nbpivot][nbr] = precision;
7724 
7725  }
7726  if(nbpivot > 0)
7727  {
7728  std::swap_ranges(S.row_begin(nbr),S.row_end(nbr),S.row_begin(nbpivot));
7729  }
7730  for(int i = nbr+1; i < M_rows; ++i)
7731  {
7732  _T u = - S[i][nbr] / S[nbr][nbr];
7733  slip::reduction(S.row_begin(i),S.row_end(i),S.row_begin(nbr),u);
7734  }
7735  nbr++;
7736  }
7737 
7738  //Calculation of X
7739  if(std::abs(S[M_rows_m1][M_cols_m1]) < precision)
7740  {
7741  ret = false;
7742  S[M.rows()-1][M.cols()-1] = precision;
7743  }
7744 
7745  bm.set_coord(0,0,M_rows_m1,M_cols_m1);
7746  slip::upper_triangular_solve(S.upper_left(bm),S.bottom_right(bm),
7747  X.begin(),X.end(),
7748  S.col_begin(M_cols),S.col_end(M_cols),
7749  precision);
7750  return ret;
7751 
7752  }
7753 
7754 /* @} */
7755 
7757  /* @{ */
7758 
7843  template<typename RandomAccessIterator1,
7844  typename RandomAccessIterator2,
7845  typename RandomAccessIterator3>
7846  void tridiagonal_lu(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
7847  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
7848  RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last,
7849  RandomAccessIterator2 L_low_first, RandomAccessIterator2 L_low_last,
7850  RandomAccessIterator3 U_up_first, RandomAccessIterator3 U_up_last,
7851  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
7852  {
7853  assert(diag_last != diag_first);
7854  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
7855  assert((diag_last - diag_first) == (low_diag_last - low_diag_first + 1));
7856  assert((L_low_last-L_low_first) == (diag_last - diag_first - 1));
7857  assert((U_up_last-U_up_first) == (diag_last - diag_first));
7858 
7859 
7860  *U_up_first = *diag_first;
7861  for(; L_low_first != L_low_last; ++L_low_first, ++U_up_first, ++diag_first, ++up_diag_first, ++low_diag_first)
7862  {
7863  if(std::abs(*U_up_first) < precision)
7864  {
7865  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
7866  }
7867  *L_low_first = *low_diag_first / *U_up_first;
7868  *(U_up_first + 1) = *(diag_first + 1) - (*L_low_first * *up_diag_first);
7869  }
7870 
7871  }
7872 
7873 /* @} */
7874 
7875 
7877  /* @{ */
7943  template<typename RandomAccessIterator1,
7944  typename RandomAccessIterator2,
7945  typename RandomAccessIterator3>
7946  void unit_lower_bidiagonal_solve(RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last,
7947  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
7948  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last)
7949  {
7950  assert(low_diag_last != low_diag_first);
7951  assert((low_diag_last - low_diag_first) == (X_last - X_first - 1));
7952  assert((low_diag_last - low_diag_first) == (B_last - B_first - 1));
7953 
7954  *X_first++ = *B_first++;
7955  for(; X_first != X_last; ++X_first, ++B_first, ++low_diag_first)
7956  {
7957  *X_first = (*B_first - (*low_diag_first * *(X_first - 1)) );
7958  }
7959  }
7960 
8005 template<typename RandomAccessIterator1,
8006  typename RandomAccessIterator2d>
8007 void unit_lower_bidiagonal_inv(RandomAccessIterator1 low_diag_first,
8008  RandomAccessIterator1 low_diag_last,
8009  RandomAccessIterator2d Ainv_up,
8010  RandomAccessIterator2d Ainv_bot)
8011 {
8012 
8013  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
8014  assert((Ainv_bot - Ainv_up)[0] == (low_diag_last - low_diag_first + 1));
8015 
8016  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8017  int n = static_cast<int>((Ainv_bot - Ainv_up)[0]);
8018  //init Ainv with the identity matrix
8019  slip::identity(Ainv_up,Ainv_bot);
8020 
8021  for(int i = 0; i < n; ++i)
8022  {
8023  slip::unit_lower_bidiagonal_solve(low_diag_first,low_diag_last,
8024  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8025  Ainv_up.col_begin(i),Ainv_up.col_end(i));
8026  }
8027 }
8028 
8099  template<typename RandomAccessIterator1, typename RandomAccessIterator2,
8100  typename RandomAccessIterator3>
8101  void lower_bidiagonal_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
8102  RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last,
8103  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
8104  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last,
8105  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
8106  {
8107  assert(diag_last != diag_first);
8108  assert((diag_last - diag_first) == (X_last - X_first));
8109  assert((diag_last - diag_first) == (B_last - B_first));
8110  assert((diag_last - diag_first) == (low_diag_last - low_diag_first + 1));
8111 
8112  if(std::abs(*diag_first) < precision)
8113  {
8114  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
8115  }
8116  *X_first++ = *B_first++ / *diag_first++;
8117  for(; X_first != X_last; ++X_first, ++B_first, ++diag_first, ++low_diag_first)
8118  {
8119  if(std::abs(*diag_first) < precision)
8120  {
8121  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
8122  }
8123  *X_first = (*B_first - (*low_diag_first * *(X_first - 1)) ) / *diag_first;
8124  }
8125  }
8126 
8127 
8180 template<typename RandomAccessIterator1,
8181  typename RandomAccessIterator2d>
8182 void lower_bidiagonal_inv(RandomAccessIterator1 diag_first,
8183  RandomAccessIterator1 diag_last,
8184  RandomAccessIterator1 low_diag_first,
8185  RandomAccessIterator1 low_diag_last,
8186  RandomAccessIterator2d Ainv_up,
8187  RandomAccessIterator2d Ainv_bot,
8188  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
8189 {
8190  assert(diag_last != diag_first);
8191  assert((diag_last - diag_first) == (low_diag_last - low_diag_first + 1));
8192  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
8193  assert((Ainv_bot - Ainv_up)[0] == (diag_last - diag_first));
8194 
8195  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8196  int n = static_cast<int>(diag_last - diag_first);
8197  //init Ainv with the identity matrix
8198  slip::identity(Ainv_up,Ainv_bot);
8199 
8200  for(int i = 0; i < n; ++i)
8201  {
8202  slip::lower_bidiagonal_solve(diag_first,diag_last,
8203  low_diag_first,low_diag_last,
8204  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8205  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8206  precision);
8207  }
8208 
8209 
8210 }
8211 
8257  template<typename Container2d1,
8258  typename Container2d2>
8259  void lower_bidiagonal_inv(const Container2d1& A,
8260  Container2d2& Ainv,
8261  typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type(1.0E-6))
8262  {
8263 
8264  assert(A.dim1() == A.dim2());
8265  assert(Ainv.dim1() == Ainv.dim2());
8266 
8267 
8268  typedef typename Container2d1::value_type value_type;
8269  typedef typename Container2d1::size_type size_type;
8270 
8271  //get the 3 main diagonals
8272  size_type n = A.dim1();
8273  slip::Array<value_type> diag(n);
8274  slip::Array<value_type> low_diag(n-1);
8275  slip::get_diagonal(A,diag.begin(),diag.end());
8276  slip::get_diagonal(A,low_diag.begin(),low_diag.end(),-1);
8277  slip::lower_bidiagonal_inv(diag.begin(),diag.end(),
8278  low_diag.begin(),low_diag.end(),
8279  Ainv.upper_left(),
8280  Ainv.bottom_right(),
8281  precision);
8282 
8283  }
8284 
8285 
8286 
8358  template<typename RandomAccessIterator1, typename RandomAccessIterator2,
8359  typename RandomAccessIterator3>
8360  void upper_bidiagonal_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
8361  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
8362  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
8363  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last,
8364  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
8365  {
8366 
8367  assert(diag_last != diag_first);
8368  assert((diag_last - diag_first) == (X_last - X_first));
8369  assert((diag_last - diag_first) == (B_last - B_first));
8370  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
8371 
8372 
8373  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8374  value_type d = *(--diag_last);
8375  if(std::abs(d) < precision)
8376  {
8377  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
8378  }
8379  *(--X_last) = *(--B_last) / d;
8380  //
8381  --X_last;
8382  --B_last;
8383  --diag_last;
8384  --up_diag_last;
8385 
8386  for(; X_last >= X_first; --X_last, --B_last, --up_diag_last, --diag_last)
8387  {
8388  if(std::abs(*diag_last) < precision)
8389  {
8390  throw std::runtime_error(slip::SINGULAR_MATRIX_ERROR);
8391  }
8392  *X_last = (*B_last - (*(X_last + 1) * *up_diag_last)) / *diag_last;
8393  }
8394 
8395  }
8396 
8449 template<typename RandomAccessIterator1,
8450  typename RandomAccessIterator2d>
8451 void upper_bidiagonal_inv(RandomAccessIterator1 diag_first,
8452  RandomAccessIterator1 diag_last,
8453  RandomAccessIterator1 up_diag_first,
8454  RandomAccessIterator1 up_diag_last,
8455  RandomAccessIterator2d Ainv_up,
8456  RandomAccessIterator2d Ainv_bot,
8457  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
8458 {
8459  assert(diag_last != diag_first);
8460  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
8461  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
8462  assert((Ainv_bot - Ainv_up)[0] == (diag_last - diag_first));
8463 
8464  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8465  int n = static_cast<int>(diag_last - diag_first);
8466  //init Ainv with the identity matrix
8467  slip::identity(Ainv_up,Ainv_bot);
8468 
8469  for(int i = 0; i < n; ++i)
8470  {
8471  slip::upper_bidiagonal_solve(diag_first,diag_last,
8472  up_diag_first,up_diag_last,
8473  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8474  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8475  precision);
8476  }
8477 
8478 
8479 }
8524  template<typename Container2d1,
8525  typename Container2d2>
8526  void upper_bidiagonal_inv(const Container2d1& A,
8527  Container2d2& Ainv,
8528  typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type(1.0E-6))
8529  {
8530 
8531  assert(A.dim1() == A.dim2());
8532  assert(Ainv.dim1() == Ainv.dim2());
8533 
8534 
8535  typedef typename Container2d1::value_type value_type;
8536  typedef typename Container2d1::size_type size_type;
8537 
8538  //get the 3 main diagonals
8539  size_type n = A.dim1();
8540  slip::Array<value_type> diag(n);
8541  slip::Array<value_type> up_diag(n-1);
8542 
8543  slip::get_diagonal(A,diag.begin(),diag.end());
8544  slip::get_diagonal(A,up_diag.begin(),up_diag.end(),1);
8545 
8546  slip::upper_bidiagonal_inv(diag.begin(),diag.end(),
8547  up_diag.begin(),up_diag.end(),
8548  Ainv.upper_left(),
8549  Ainv.bottom_right(),
8550  precision);
8551 
8552  }
8553 
8618  template<typename RandomAccessIterator1, typename RandomAccessIterator2,
8619  typename RandomAccessIterator3>
8620  void unit_upper_bidiagonal_solve(RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
8621  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
8622  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last
8623  )
8624  {
8625 
8626  assert(X_last != X_first);
8627  assert((X_last - X_first) == (B_last - B_first));
8628  assert((X_last - X_first) == (up_diag_last - up_diag_first + 1));
8629 
8630 
8631  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8632  *(--X_last) = *(--B_last);
8633  //
8634  --X_last;
8635  --B_last;
8636  --up_diag_last;
8637 
8638  for(; X_last >= X_first; --X_last, --B_last, --up_diag_last)
8639  {
8640  *X_last = (*B_last - (*(X_last + 1) * *up_diag_last));
8641  }
8642  }
8643 
8644 
8688 template<typename RandomAccessIterator1,
8689  typename RandomAccessIterator2d>
8690 void unit_upper_bidiagonal_inv(RandomAccessIterator1 up_diag_first,
8691  RandomAccessIterator1 up_diag_last,
8692  RandomAccessIterator2d Ainv_up,
8693  RandomAccessIterator2d Ainv_bot)
8694 {
8695  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
8696  assert((Ainv_bot - Ainv_up)[0] == (up_diag_last - up_diag_first + 1));
8697 
8698  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8699  int n = static_cast<int>((Ainv_bot - Ainv_up)[0]);
8700  //init Ainv with the identity matrix
8701  slip::identity(Ainv_up,Ainv_bot);
8702 
8703  for(int i = 0; i < n; ++i)
8704  {
8705  slip::unit_upper_bidiagonal_solve(up_diag_first,up_diag_last,
8706  Ainv_up.col_begin(i),Ainv_up.col_end(i),
8707  Ainv_up.col_begin(i),Ainv_up.col_end(i));
8708  }
8709 
8710 
8711 }
8712 /* @} */
8713 
8715  /* @{ */
8792  template<typename RandomAccessIterator1,
8793  typename RandomAccessIterator2,
8794  typename RandomAccessIterator3>
8795  void thomas_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
8796  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
8797  RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last,
8798  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
8799  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last,
8800  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
8801  {
8802  assert(diag_last != diag_first);
8803  assert((diag_last - diag_first) == (X_last - X_first));
8804  assert((diag_last - diag_first) == (B_last - B_first));
8805  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
8806  assert((low_diag_last - low_diag_first) == (up_diag_last - up_diag_first));
8807 
8808  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference;
8809  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
8810  _Difference n = (diag_last - diag_first);
8811 
8812  //tridiagonal LU decomposition
8813  slip::Array<value_type> M(n);
8814  slip::Array<value_type> L(n-1);
8815  tridiagonal_lu(diag_first,diag_last,
8816  up_diag_first,up_diag_last,
8817  low_diag_first,low_diag_last,
8818  L.begin(),L.end(),
8819  M.begin(),M.end(),precision);
8820 
8821  //First step: resolution of the lower bidiagonal system LY = B
8822  //L main diagonal are ones
8823  //L lower diagonal is L vector
8824  slip::Array<value_type> Y(n);
8825  unit_lower_bidiagonal_solve(L.begin(),L.end(),Y.begin(),Y.end(),B_first,B_last);
8826  //Second step: resolution of the upper bidiagonal system UX = Y
8827  //U main diagonal is M vector
8828  //U upper diagonal is up_diag range
8829  upper_bidiagonal_solve(M.begin(),M.end(),up_diag_first,up_diag_last,X_first,X_last,Y.begin(),Y.end(),precision);
8830 
8831  }
8832 
8833 
8908  template<typename Container2d,
8909  typename Vector1,
8910  typename Vector2>
8911  void thomas_solve(const Container2d& A,
8912  Vector1& X,
8913  const Vector2& B,
8914  typename slip::lin_alg_traits<typename Container2d::value_type>::value_type precision = typename slip::lin_alg_traits<typename Container2d::value_type>::value_type(1.0E-6))
8915  {
8916 
8917  assert(A.dim1() == A.dim2());
8918  assert(X.size() == A.dim1());
8919  assert(B.size() == A.dim1());
8920 
8921  typedef typename Container2d::value_type value_type;
8922  typedef typename Container2d::size_type size_type;
8923 
8924  //get the 3 main diagonals
8925  size_type n = A.dim1();
8926  slip::Array<value_type> diag(n);
8927  slip::Array<value_type> up_diag(n-1);
8928  slip::Array<value_type> low_diag(n-1);
8929  slip::get_diagonal(A,diag.begin(),diag.end());
8930  slip::get_diagonal(A,up_diag.begin(),up_diag.end(),1);
8931  slip::get_diagonal(A,low_diag.begin(),low_diag.end(),-1);
8932  //tridiagonal solve
8933  thomas_solve(diag.begin(),diag.end(),
8934  up_diag.begin(),up_diag.end(),
8935  low_diag.begin(),low_diag.end(),
8936  X.begin(),X.end(),
8937  B.begin(),B.end(),precision);
8938  }
8939 
8998 template<typename RandomAccessIterator1,
8999  typename RandomAccessIterator2d>
9000 void thomas_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
9001  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
9002  RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last,
9003  RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot,
9004  typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type precision = typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type(1.0E-6))
9005 {
9006  assert(diag_last != diag_first);
9007  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
9008  assert((low_diag_last - low_diag_first) == (up_diag_last - up_diag_first));
9009  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
9010  assert((Ainv_bot - Ainv_up)[0] == (diag_last - diag_first));
9011  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference;
9012  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
9013  _Difference n = (diag_last - diag_first);
9014  //init Ainv with the identity matrix
9015  slip::identity(Ainv_up,Ainv_bot);
9016 
9017  for(_Difference i = 0; i < n; ++i)
9018  {
9019  slip::thomas_solve(diag_first,diag_last,
9020  up_diag_first,up_diag_last,
9021  low_diag_first,low_diag_last,
9022  Ainv_up.col_begin(i),Ainv_up.col_end(i),
9023  Ainv_up.col_begin(i),Ainv_up.col_end(i),
9024  precision);
9025  }
9026 
9027 
9028 }
9076  template<typename Container2d1,
9077  typename Container2d2>
9078  void thomas_inv(const Container2d1& A,
9079  Container2d2& Ainv,
9080  typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Container2d1::value_type>::value_type(1.0E-6))
9081  {
9082 
9083  assert(A.dim1() == A.dim2());
9084  assert(Ainv.dim1() == Ainv.dim2());
9085 
9086 
9087  typedef typename Container2d1::value_type value_type;
9088  typedef typename Container2d1::size_type size_type;
9089 
9090  //get the 3 main diagonals
9091  size_type n = A.dim1();
9092  slip::Array<value_type> diag(n);
9093  slip::Array<value_type> up_diag(n-1);
9094  slip::Array<value_type> low_diag(n-1);
9095  slip::get_diagonal(A,diag.begin(),diag.end());
9096  slip::get_diagonal(A,up_diag.begin(),up_diag.end(),1);
9097  slip::get_diagonal(A,low_diag.begin(),low_diag.end(),-1);
9098  slip::thomas_inv(diag.begin(),diag.end(),
9099  up_diag.begin(),up_diag.end(),
9100  low_diag.begin(),low_diag.end(),
9101  Ainv.upper_left(),
9102  Ainv.bottom_right(),
9103  precision);
9104 
9105  }
9106 /* @} */
9107 
9109  /* @{ */
9110 
9200 template<typename RandomAccessIterator1,
9201  typename RandomAccessIterator2,
9202  typename RandomAccessIterator3>
9203 void tridiagonal_cholesky(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
9204  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
9205 
9206  RandomAccessIterator2 R_diag_first, RandomAccessIterator2 R_diag_last,
9207  RandomAccessIterator3 R_low_first, RandomAccessIterator3 R_low_last)
9208  {
9209  assert(diag_last != diag_first);
9210  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
9211  assert((R_diag_last-R_diag_first) == (diag_last - diag_first));
9212  assert((R_diag_last-R_diag_first) == (R_low_last - R_low_first + 1));
9213 
9214  typedef typename slip::lin_alg_traits<typename std::iterator_traits<RandomAccessIterator1>::value_type>::value_type real_type;
9215  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
9216  if(std::real(*diag_first) <= real_type(0))
9217  {
9218  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
9219  }
9220  *R_diag_first = std::sqrt(*diag_first);
9221 
9222  for(; R_low_first != R_low_last; ++R_low_first, ++up_diag_first, ++diag_first, ++R_diag_first)
9223  {
9224 
9225 
9226  *R_low_first = slip::conj(*up_diag_first / *R_diag_first);
9227  value_type proj = *(diag_first + 1) - slip::conj(*R_low_first) * *R_low_first;
9228  if(std::real(proj) <= real_type(0))
9229  {
9230  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
9231  }
9232  *(R_diag_first + 1) = std::sqrt(proj);
9233  }
9234 
9235  }
9318 template<typename RandomAccessIterator1,
9319  typename RandomAccessIterator2,
9320  typename RandomAccessIterator3>
9321  void tridiagonal_cholesky_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
9322  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
9323  RandomAccessIterator2 X_first, RandomAccessIterator2 X_last,
9324  RandomAccessIterator3 B_first, RandomAccessIterator3 B_last)
9325  {
9326  assert((X_last - X_first) == (B_last - B_first));
9327  assert((diag_last - diag_first) == (X_last - X_first));
9328  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
9329  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference;
9330  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
9331  _Difference n = (diag_last - diag_first);
9332 
9333  //tridiagonal LU decomposition
9334  slip::Array<value_type> R_diag(n);
9335  slip::Array<value_type> R_low(n-1);
9336  slip::tridiagonal_cholesky(diag_first,diag_last,
9337  up_diag_first,up_diag_last,
9338  R_diag.begin(),R_diag.end(),
9339  R_low.begin(),R_low.end());
9340 
9341  //First step: resolution of the lower bidiagonal system LY = B
9342  //L main diagonal are ones
9343  //L lower diagonal is L vector
9344  slip::Array<value_type> Y(n);
9345  slip::lower_bidiagonal_solve(R_diag.begin(),R_diag.end(),
9346  R_low.begin(),R_low.end(),
9347  Y.begin(),Y.end(),
9348  B_first,B_last);
9349  //Second step: resolution of the upper bidiagonal system UX = Y
9350  //U main diagonal is M vector
9351  //U upper diagonal is up_diag range
9352  slip::apply(R_low.begin(),R_low.end(),R_low.begin(),slip::conj);
9353  slip::upper_bidiagonal_solve(R_diag.begin(),R_diag.end(),
9354  R_low.begin(),R_low.end(),
9355  X_first,X_last,
9356  Y.begin(),Y.end());
9357 
9358  }
9359 
9434 template<typename Container2d,
9435  typename Vector1,
9436  typename Vector2>
9437  void tridiagonal_cholesky_solve(const Container2d& A,
9438  Vector1& X,
9439  const Vector2& B)
9440  {
9441 
9442  assert(A.dim1() == A.dim2());
9443  assert(X.size() == A.dim1());
9444  assert(B.size() == A.dim1());
9445 
9446  typedef typename Container2d::value_type value_type;
9447  typedef typename Container2d::size_type size_type;
9448 
9449  //get the 2 main diagonals
9450  size_type n = A.dim1();
9451  slip::Array<value_type> diag(n);
9452  slip::Array<value_type> up_diag(n-1);
9453  slip::get_diagonal(A,diag.begin(),diag.end());
9454  slip::get_diagonal(A,up_diag.begin(),up_diag.end(),1);
9455  //tridiagonal solve
9456  slip::tridiagonal_cholesky_solve(diag.begin(),diag.end(),
9457  up_diag.begin(),up_diag.end(),
9458  X.begin(),X.end(),
9459  B.begin(),B.end());
9460  }
9461 
9514  template<typename RandomAccessIterator1,
9515  typename RandomAccessIterator2d>
9516 void tridiagonal_cholesky_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last,
9517  RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last,
9518  RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot)
9519 {
9520  assert(diag_last != diag_first);
9521  assert((diag_last - diag_first) == (up_diag_last - up_diag_first + 1));
9522  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
9523  assert((Ainv_bot - Ainv_up)[0] == (diag_last - diag_first));
9524  typedef typename std::iterator_traits<RandomAccessIterator1>::difference_type _Difference;
9525  typedef typename std::iterator_traits<RandomAccessIterator1>::value_type value_type;
9526  _Difference n = (diag_last - diag_first);
9527  //init Ainv with the identity matrix
9528  slip::identity(Ainv_up,Ainv_bot);
9529  for(_Difference i = 0; i < n; ++i)
9530  {
9531  slip::tridiagonal_cholesky_solve(diag_first,diag_last,
9532  up_diag_first,up_diag_last,
9533  Ainv_up.col_begin(i),Ainv_up.col_end(i),
9534  Ainv_up.col_begin(i),Ainv_up.col_end(i));
9535  }
9536 
9537 
9538 }
9539 
9587 template<typename Container2d1,
9588  typename Container2d2>
9589 void tridiagonal_cholesky_inv(const Container2d1& A,
9590  Container2d2& Ainv)
9591  {
9592 
9593  assert(A.dim1() == A.dim2());
9594  assert(Ainv.dim1() == Ainv.dim2());
9595 
9596 
9597  typedef typename Container2d1::value_type value_type;
9598  typedef typename Container2d1::size_type size_type;
9599 
9600  //get the 3 main diagonals
9601  size_type n = A.dim1();
9602  slip::Array<value_type> diag(n);
9603  slip::Array<value_type> up_diag(n-1);
9604 
9605  slip::get_diagonal(A,diag.begin(),diag.end());
9606  slip::get_diagonal(A,up_diag.begin(),up_diag.end(),1);
9607 
9608  slip::tridiagonal_cholesky_inv(diag.begin(),diag.end(),
9609  up_diag.begin(),up_diag.end(),
9610  Ainv.upper_left(),
9611  Ainv.bottom_right());
9612 
9613  }
9614 
9615 /* @} */
9616 
9618  /* @{ */
9619 
9632  template<typename Container2d>
9633  void LDLT_decomposition(Container2d& A)
9634  {
9635  typedef typename Container2d::value_type value_type;
9636  typedef typename Container2d::size_type size_type;
9637  int n = int(A.dim1());
9638  slip::Array<value_type> R(n);
9639  value_type proj = value_type();
9640  for(int i = 0; i < n; ++i)
9641  {
9642  //computes R
9643  for(int j = 0; j < i; ++j)
9644  {
9645  R[j] = A[i][j] * A[j][j];
9646  }
9647 
9648  proj = std::inner_product(A.row_begin(i),A.row_begin(i)+i,
9649  R.begin(),
9650  value_type());
9651  R[i] = A[i][i] - proj;
9652  //store D
9653  A[i][i] = R[i];
9654  //computes L
9655  for(int k = (i+1); k < n; ++k)
9656  {
9657  proj = std::inner_product(A.row_begin(k),A.row_begin(k)+i,
9658  R.begin(),
9659  value_type());
9660  A[k][i] = (A[k][i] - proj) / R[i];
9661  }
9662  }
9663 
9664  }
9665 
9666 
9679  template<typename Matrix1,
9680  typename Matrix2,
9681  typename Vector1,
9682  typename Matrix3>
9683  void LDLT_decomposition(const Matrix1& A,
9684  Matrix2& L,
9685  Vector1& D,
9686  Matrix3& LT)
9687  {
9688  typedef typename Matrix1::value_type value_type;
9689  typedef typename Matrix1::size_type size_type;
9690  int n = int(A.dim1());
9691  slip::Array<value_type> R(n);
9692  value_type proj = value_type();
9693  for(int i = 0; i < n; ++i)
9694  {
9695  //computes R
9696  for(int j = 0; j < i; ++j)
9697  {
9698  R[j] = L[i][j] * D[j];
9699  }
9700 
9701  proj = std::inner_product(L.row_begin(i),L.row_begin(i)+i,
9702  R.begin(),
9703  value_type());
9704  R[i] = A[i][i] - proj;
9705  D[i] = R[i];
9706  L[i][i] = value_type(1);
9707  LT[i][i] = value_type(1);
9708  //store D
9709  //A[i][i] = R[i];
9710  //computes L
9711  for(int k = (i+1); k < n; ++k)
9712  {
9713  proj = std::inner_product(L.row_begin(k),L.row_begin(k)+i,
9714  R.begin(),
9715  value_type());
9716  L[k][i] = (A[k][i] - proj) / R[i];
9717  LT[i][k] = L[k][i];
9718  }
9719  }
9720 
9721  }
9722 
9735  template<typename Matrix1,
9736  typename Vector1,
9737  typename Vector2>
9738  void LDLT_solve(const Matrix1& A,
9739  Vector1& X,
9740  const Vector2& B)
9741  {
9742  typedef typename Matrix1::value_type value_type;
9743  typedef typename Matrix1::size_type size_type;
9744  size_type n = A.dim1();
9745  //decomposition
9746  slip::Array2d<value_type> L(n,n);
9747  slip::Array<value_type> D(n);
9748  slip::Array2d<value_type> LT(n,n);
9749  slip::LDLT_decomposition(A,L,D,LT);
9750 
9751  //resolution of LY = B
9752  slip::Array<value_type> Y(n);
9753  slip::unit_lower_triangular_solve(L,Y,B);
9754  //resolution of DZ = Y
9755  slip::Array<value_type> Z(n);
9757  slip::diagonal_solve(D.begin(),D.end(),
9758  Z.begin(),Z.end(),
9759  Y.begin(),Y.end(),1E-06);
9760 
9761  //resolution of LTX = Z
9762  slip::unit_upper_triangular_solve(LT,X,Z);
9763  }
9764 /* @} */
9765 
9767  /* @{ */
9818  template <typename MatrixIterator1,
9819  typename MatrixIterator2,
9820  typename MatrixIterator3>
9821  inline
9822  void cholesky_cplx(MatrixIterator1 A_up,
9823  MatrixIterator1 A_bot,
9824  MatrixIterator2 L_up,
9825  MatrixIterator2 L_bot,
9826  MatrixIterator3 LH_up,
9827  MatrixIterator3 LH_bot)
9828  {
9829  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
9830  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
9831  assert((LH_bot - LH_up)[0] == (LH_bot - LH_up)[1]);
9832  assert((A_bot - A_up)[0] == (L_bot - L_up)[0]);
9833  assert((A_bot - A_up)[0] == (LH_bot - LH_up)[0]);
9834  assert((L_bot - L_up)[0] == (LH_bot - LH_up)[0]);
9835 
9836  typedef typename MatrixIterator1::value_type value_type;
9837  typedef typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type real_type;
9838  std::size_t A_rows = std::size_t((A_bot - A_up)[0]);
9839 
9840  slip::DPoint2d<int> d(0,0);
9841  slip::DPoint2d<int> d2(0,0);
9842  slip::DPoint2d<int> d3(0,0);
9843  value_type proj = value_type();
9844  for(std::size_t i = 0; i < A_rows; ++i)
9845  {
9846  d.set_coord(i,i);
9847  proj = slip::hermitian_inner_product(L_up.row_begin(i),
9848  L_up.row_begin(i)+i,
9849  L_up.row_begin(i),
9850  value_type());
9851  value_type Aup_m_proj = A_up[d] - proj;
9852  if(std::real(Aup_m_proj) <= real_type())
9853  {
9854  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
9855  }
9856  value_type Lii = std::sqrt(Aup_m_proj);
9857  L_up[d] = Lii;
9858  LH_up[d] = Lii;
9859  for(std::size_t j = (i+1); j < A_rows; ++j)
9860  {
9861  proj = slip::hermitian_inner_product(L_up.row_begin(j),
9862  L_up.row_begin(j)+i,
9863  L_up.row_begin(i),
9864  value_type());
9865  d2.set_coord(j,i);
9866  d3.set_coord(i,j);
9867  L_up[d2] = (A_up[d2] - slip::conj(proj))/Lii;
9868  L_up[d3] = value_type();
9869  LH_up[d3] = slip::conj(L_up[d2]);
9870  LH_up[d2] = value_type();
9871  }
9872  }
9873 
9874 
9875  }
9876 
9877 
9911  template <typename MatrixIterator1>
9912  inline
9913  void cholesky_cplx(MatrixIterator1 A_up,
9914  MatrixIterator1 A_bot)
9915 
9916  {
9917  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
9918 
9919  typedef typename MatrixIterator1::value_type value_type;
9920  typedef typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type real_type;
9921  std::size_t A_rows = std::size_t((A_bot - A_up)[0]);
9922 
9923  slip::DPoint2d<int> d(0,0);
9924  slip::DPoint2d<int> d2(0,0);
9925  slip::DPoint2d<int> d3(0,0);
9926  value_type proj = value_type();
9927  for(std::size_t i = 0; i < A_rows; ++i)
9928  {
9929  d.set_coord(i,i);
9930  proj = slip::hermitian_inner_product(A_up.row_begin(i),
9931  A_up.row_begin(i)+i,
9932  A_up.row_begin(i),
9933  value_type());
9934  value_type Aup_m_proj = A_up[d] - proj;
9935  if(std::real(Aup_m_proj) <= real_type())
9936  {
9937  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
9938  }
9939  value_type Aii = std::sqrt(Aup_m_proj);
9940 
9941  A_up[d] = Aii;
9942  for(std::size_t j = (i+1); j < A_rows; ++j)
9943  {
9944  proj = slip::hermitian_inner_product(A_up.row_begin(j),
9945  A_up.row_begin(j)+i,
9946  A_up.row_begin(i),
9947  value_type());
9948  d2.set_coord(j,i);
9949  d3.set_coord(i,j);
9950  A_up[d2] = (A_up[d2] - slip::conj(proj))/Aii;
9951  A_up[d3] = slip::conj(A_up[d2]);
9952 
9953  }
9954  }
9955 
9956  }
9957 
10004  template <typename MatrixIterator1,
10005  typename MatrixIterator2,
10006  typename MatrixIterator3>
10007  inline
10008  void cholesky_real(MatrixIterator1 A_up,
10009  MatrixIterator1 A_bot,
10010  MatrixIterator2 L_up,
10011  MatrixIterator2 L_bot,
10012  MatrixIterator3 LH_up,
10013  MatrixIterator3 LH_bot)
10014  {
10015  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
10016  assert((L_bot - L_up)[0] == (L_bot - L_up)[1]);
10017  assert((LH_bot - LH_up)[0] == (LH_bot - LH_up)[1]);
10018  assert((A_bot - A_up)[0] == (L_bot - L_up)[0]);
10019  assert((A_bot - A_up)[0] == (LH_bot - LH_up)[0]);
10020  assert((L_bot - L_up)[0] == (LH_bot - LH_up)[0]);
10021 
10022  typedef typename MatrixIterator1::value_type value_type;
10023  typedef typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type real_type;
10024  std::size_t A_rows = std::size_t((A_bot - A_up)[0]);
10025 
10026  slip::DPoint2d<int> d(0,0);
10027  slip::DPoint2d<int> d2(0,0);
10028  slip::DPoint2d<int> d3(0,0);
10029  value_type proj = value_type();
10030  for(std::size_t i = 0; i < A_rows; ++i)
10031  {
10032  d.set_coord(i,i);
10033  proj = slip::inner_product(L_up.row_begin(i),
10034  L_up.row_begin(i)+i,
10035  L_up.row_begin(i),
10036  value_type());
10037  value_type Aup_m_proj = A_up[d] - proj;
10038  if(Aup_m_proj <= real_type())
10039  {
10040  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
10041  }
10042  value_type Lii = std::sqrt(Aup_m_proj);
10043 
10044  L_up[d] = Lii;
10045  LH_up[d] = Lii;
10046  for(std::size_t j = (i+1); j < A_rows; ++j)
10047  {
10048  proj = slip::inner_product(L_up.row_begin(j),
10049  L_up.row_begin(j)+i,
10050  L_up.row_begin(i),
10051  value_type());
10052  d2.set_coord(j,i);
10053  d3.set_coord(i,j);
10054  L_up[d2] = (A_up[d2] - proj)/Lii;
10055  L_up[d3] = value_type();
10056  LH_up[d3] = L_up[d2];
10057  LH_up[d2] = value_type();
10058  }
10059  }
10060 
10061  }
10091  template <typename MatrixIterator1>
10092  inline
10093  void cholesky_real(MatrixIterator1 A_up,
10094  MatrixIterator1 A_bot)
10095  {
10096  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
10097 
10098  typedef typename MatrixIterator1::value_type value_type;
10099  typedef typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type real_type;
10100  std::size_t A_rows = std::size_t((A_bot - A_up)[0]);
10101 
10102  //init L and LH
10103  slip::DPoint2d<int> d(0,0);
10104  slip::DPoint2d<int> d2(0,0);
10105  slip::DPoint2d<int> d3(0,0);
10106  value_type proj = value_type();
10107  for(std::size_t i = 0; i < A_rows; ++i)
10108  {
10109  d.set_coord(i,i);
10110  proj = slip::inner_product(A_up.row_begin(i),
10111  A_up.row_begin(i)+i,
10112  A_up.row_begin(i),
10113  value_type());
10114  value_type Aup_m_proj = A_up[d] - proj;
10115  if(std::real(Aup_m_proj) <= real_type(0))
10116  {
10117  throw std::runtime_error(slip::POSITIVE_DEFINITE_MATRIX_ERROR);
10118  }
10119  value_type Aii = std::sqrt(Aup_m_proj);
10120 
10121  A_up[d] = Aii;
10122  for(std::size_t j = (i+1); j < A_rows; ++j)
10123  {
10124  proj = slip::inner_product(A_up.row_begin(j),
10125  A_up.row_begin(j)+i,
10126  A_up.row_begin(i),
10127  value_type());
10128  d2.set_coord(j,i);
10129  d3.set_coord(i,j);
10130  A_up[d2] = (A_up[d2] - proj)/Aii;
10131  A_up[d3] = A_up[d2];
10132  }
10133  }
10134 
10135 
10136  }
10137 
10138 
10176  template <class Matrix1,
10177  class Matrix2,
10178  class Matrix3>
10179  inline
10180  void cholesky_real(const Matrix1 & A,
10181  Matrix2 & L,
10182  Matrix3& LT)
10183  {
10184  slip::cholesky_real(A.upper_left(),A.bottom_right(),
10185  L.upper_left(),L.bottom_right(),
10186  LT.upper_left(),LT.bottom_right());
10187  }
10188 
10189 
10231  template <class Matrix1,
10232  class Matrix2,
10233  class Matrix3>
10234  inline
10235  void cholesky_cplx(const Matrix1 & A,
10236  Matrix2 & L,
10237  Matrix3& LT)
10238  {
10239  slip::cholesky_cplx(A.upper_left(),A.bottom_right(),
10240  L.upper_left(),L.bottom_right(),
10241  LT.upper_left(),LT.bottom_right());
10242  }
10243 
10244 
10286  template <class Matrix1,
10287  class Matrix2,
10288  class Matrix3>
10289  inline
10290  void cholesky(const Matrix1 & A,
10291  Matrix2 & L,
10292  Matrix3& LT)
10293  {
10294  if(slip::is_complex(typename Matrix1::value_type()))
10295  {
10296  slip::cholesky_cplx(A.upper_left(),A.bottom_right(),
10297  L.upper_left(),L.bottom_right(),
10298  LT.upper_left(),LT.bottom_right());
10299  }
10300  else
10301  {
10302  slip::cholesky_real(A.upper_left(),A.bottom_right(),
10303  L.upper_left(),L.bottom_right(),
10304  LT.upper_left(),LT.bottom_right());
10305  }
10306  }
10307 
10346  template <class Matrix1>
10347  inline
10348  void cholesky(Matrix1 & A)
10349 
10350  {
10351  if(slip::is_complex(typename Matrix1::value_type()))
10352  {
10353  slip::cholesky_cplx(A.upper_left(),A.bottom_right());
10354  }
10355  else
10356  {
10357  slip::cholesky_real(A.upper_left(),A.bottom_right());
10358  }
10359  }
10360 
10408  template <typename MatrixIterator,
10409  typename RandomAccessIterator1,
10410  typename RandomAccessIterator2>
10411  inline
10412  void cholesky_solve(MatrixIterator A_up,
10413  MatrixIterator A_bot,
10414  RandomAccessIterator1 X_first,
10415  RandomAccessIterator1 X_last,
10416  RandomAccessIterator2 B_first,
10417  RandomAccessIterator2 B_last,
10418  typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type precision = typename slip::lin_alg_traits<typename MatrixIterator::value_type>::value_type(1.0E-6))
10419  {
10420  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
10421  assert((A_bot - A_up)[1] == (X_last - X_first));
10422  assert((X_last - X_first) == (B_last - B_first));
10423 
10424  typedef typename std::iterator_traits<MatrixIterator>::value_type value_type;
10425  typedef typename MatrixIterator::size_type size_type;
10426  size_type n = size_type((A_bot - A_up)[0]);
10427 
10428  //decomposition
10429  slip::Array2d<value_type> LLT(n,n);
10430  std::copy(A_up,A_bot,LLT.begin());
10431  slip::cholesky(LLT);
10432  slip::Array<value_type> Y(n);
10433  //resolution of LY = B
10434  slip::lower_triangular_solve(LLT.upper_left(),LLT.bottom_right(),
10435  Y.begin(),Y.end(),
10436  B_first,B_last,precision);
10437  //resolution of LTX = Y
10438  slip::upper_triangular_solve(LLT.upper_left(),LLT.bottom_right(),
10439  X_first,X_last,
10440  Y.begin(),Y.end(),precision);
10441 
10442  }
10443 
10444 
10486  template <class Matrix1,
10487  class Vector1,
10488  class Vector2>
10489  inline
10490  void cholesky_solve(const Matrix1 & A,
10491  Vector1 & X,
10492  const Vector2& B,
10493  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type precision =
10494  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type(1.0E-6))
10495  {
10496  slip::cholesky_solve(A.upper_left(),A.bottom_right(),
10497  X.begin(),X.end(),
10498  B.begin(),B.end(),
10499  precision);
10500 
10501  }
10502 
10503 
10547  template <typename MatrixIterator1,
10548  typename MatrixIterator2>
10549  inline
10550  void cholesky_inv(MatrixIterator1 A_up,
10551  MatrixIterator1 A_bot,
10552  MatrixIterator2 Ainv_up,
10553  MatrixIterator2 Ainv_bot,
10554  typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type precision = typename slip::lin_alg_traits<typename MatrixIterator1::value_type>::value_type(1.0E-6))
10555  {
10556  assert((A_bot - A_up)[0] == (A_bot - A_up)[1]);
10557  assert((Ainv_bot - Ainv_up)[0] == (Ainv_bot - Ainv_up)[1]);
10558  assert((A_bot - A_up)[0] == (Ainv_bot - Ainv_up)[0]);
10559  assert((A_bot - A_up)[1] == (Ainv_bot - Ainv_up)[1]);
10560 
10561 
10562  typedef typename MatrixIterator1::value_type value_type;
10563  std::size_t n = std::size_t((A_bot - A_up)[0]);
10564  //init Ainv with the identity matrix
10565  slip::identity(Ainv_up,Ainv_bot);
10566 
10567  slip::Array2d<value_type> LLT(n,n);
10568  std::copy(A_up,A_bot,LLT.upper_left());
10569  //cholesky decomposition
10570  if(slip::is_complex(value_type()))
10571  {
10572  slip::cholesky_cplx(LLT.upper_left(),LLT.bottom_right());
10573  }
10574  else
10575  {
10576  slip::cholesky_real(LLT.upper_left(),LLT.bottom_right());
10577  }
10578 
10579  //solve the n systems
10580  slip::Array<value_type> Y(n);
10581  for(std::size_t k = 0; k < n; ++k)
10582  {
10583  //resolution of LY = B
10584  slip::lower_triangular_solve(LLT.upper_left(),LLT.bottom_right(),
10585  Y.begin(),Y.end(),
10586  Ainv_up.col_begin(k),Ainv_up.col_end(k),
10587  precision);
10588  //resolution of LTX = Y
10589  slip::upper_triangular_solve(LLT.upper_left(),LLT.bottom_right(),
10590  Ainv_up.col_begin(k),Ainv_up.col_end(k),
10591  Y.begin(),Y.end(),precision);
10592  }
10593 
10594  }
10595 
10634 template <typename Matrix1,
10635  typename Matrix2>
10636  inline
10637  void cholesky_inv(Matrix1 A,
10638  Matrix2 Ainv,
10639  typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type precision = typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type(1.0E-6))
10640  {
10641  slip::cholesky_inv(A.upper_left(),A.bottom_right(),
10642  Ainv.upper_left(),Ainv.bottom_right(),
10643  precision);
10644  }
10645  /* @} */
10646 
10647 
10649 /* @{ */
10662 template <typename Real >
10663  inline
10664  void rotgen(Real& a,
10665  Real& b,
10666  Real& cos,
10667  Real& sin)
10668  {
10669  if(b == Real(0.0))
10670  {
10671  cos = Real(1.0);
10672  sin = Real(0.0);
10673  return;
10674  }
10675  if(a == Real(0.0))
10676  {
10677  cos = Real(0.0);
10678  sin = Real(1.0);
10679  a = b;
10680  b = Real(0.0);
10681  return;
10682  }
10683  Real mu = a/(std::abs(a));
10684  Real tau = std::abs(a)+std::abs(b);
10685  Real nu = tau*(std::sqrt(std::abs(a/tau)*std::abs(a/tau)+std::abs(b/tau)*std::abs(b/tau)));
10686  cos = std::abs(a)/nu;
10687  sin = mu * ((slip::conj(b))/nu);
10688  a = nu * mu;
10689  b = Real(0.0);
10690  }
10691 
10718  template <typename Real>
10719  inline
10720  void givens_sin_cos(const Real& xi, const Real& xk,
10721  Real& sin, Real& cos)
10722  {
10723  Real n = std::sqrt(xi*xi + xk * xk);
10724  cos = xi /n;
10725  sin = -xk/n;
10726  }
10727 
10728 
10759  template <typename Complex>
10760  inline
10761  void complex_givens_sin_cos(const Complex& xi, const Complex& xk,
10762  Complex& sin, typename Complex::value_type& cos)
10763  {
10764  typedef typename Complex::value_type Real;
10765  Real epsilon = slip::sign(xi.real());
10766  Real n = std::sqrt(std::norm(xi) + std::norm(xk));
10767  cos = std::abs(xi) / (epsilon * n);
10768  sin = - xk / (epsilon * (slip::complex_sign<Complex>()(xi) * n));
10769  }
10770 
10771 
10772 
10796  template <typename Matrix, typename SizeType, typename Real>
10797  inline
10798  void left_givens(Matrix & M,
10799  const SizeType& row1,
10800  const SizeType& row2,
10801  const Real& sinus,
10802  const Real& cosinus,
10803  const SizeType& col1,
10804  const SizeType& col2)
10805  {
10806  typedef typename Matrix::value_type value_type;
10807  for(SizeType j = col1; j <= col2; ++j)
10808  {
10809  value_type t1 = M(row1,j);
10810  value_type t2 = M(row2,j);
10811  M(row1,j) = cosinus * t1 + sinus * t2;
10812  M(row2,j) = - sinus * t1 + cosinus * t2;
10813  }
10814  }
10815 
10839  template <typename Matrix, typename SizeType, typename Real>
10840  inline
10841  void right_givens(Matrix & M,
10842  const SizeType& col1,
10843  const SizeType& col2,
10844  const Real& sinus,
10845  const Real& cosinus,
10846  const SizeType& row1,
10847  const SizeType& row2)
10848  {
10849  typedef typename Matrix::value_type value_type;
10850  for(SizeType i = row1; i <= row2; ++i)
10851  {
10852  value_type t1 = M(i,col1);
10853  value_type t2 = M(i,col2);
10854  M(i,col1) = cosinus * t1 + sinus * t2;
10855  M(i,col2) = - sinus * t1 + cosinus * t2;
10856  }
10857  }
10858 
10859 
10883  template <typename Matrix, typename SizeType, typename Complex>
10884  inline
10885  void complex_left_givens(Matrix & M,
10886  const SizeType& row1,
10887  const SizeType& row2,
10888  const Complex& sinus,
10889  const typename Complex::value_type& cosinus,
10890  const SizeType& col1,
10891  const SizeType& col2)
10892  {
10893  for(SizeType j = col1; j <= col2; ++j)
10894  {
10895  Complex t1 = M(row1,j);
10896  Complex t2 = M(row2,j);
10897  M(row1,j) = cosinus * t1 + slip::conj(sinus) * t2;
10898  M(row2,j) = -sinus * t1 + cosinus * t2;
10899  }
10900  }
10901 
10902 
10926  template <typename Matrix, typename SizeType, typename Complex>
10927  inline
10928  void complex_right_givens(Matrix & M,
10929  const SizeType& col1,
10930  const SizeType& col2,
10931  const Complex& sinus,
10932  const typename Complex::value_type& cosinus,
10933  const SizeType& row1,
10934  const SizeType& row2)
10935  {
10936  for(SizeType i = row1; i <= row2; ++i)
10937  {
10938  Complex t1 = M(i,col1);
10939  Complex t2 = M(i,col2);
10940  M(i,col1) = cosinus * t1 + sinus * t2;
10941  M(i,col2) = -slip::conj(sinus) * t1 + cosinus * t2;
10942  }
10943  }
10944 
10945 
10947 /* @{ */
10948 
10979 template <typename VectorIterator1,typename VectorIterator2>
10980 inline
10981 void housegen(VectorIterator1 a_begin, VectorIterator1 a_end,
10982  VectorIterator2 u_begin,VectorIterator2 u_end,
10983  typename std::iterator_traits<VectorIterator1>::value_type &nu)
10984 
10985  {
10986 
10987  assert((a_end - a_begin) == (u_end - u_begin));
10988  typedef typename std::iterator_traits<VectorIterator1>::value_type value_type;
10989  typedef typename slip::lin_alg_traits<typename std::iterator_traits<VectorIterator1>::value_type>::value_type norm_type;
10990  value_type rho = value_type(0);
10991 
10992  std::copy(a_begin,a_end,u_begin);
10993  nu = std::sqrt(slip::L22_norm_cplx(a_begin,a_end));
10994  if(std::abs(nu) == value_type(0.0))
10995  {
10996  *(u_begin) = value_type(std::sqrt(2.0));
10997  return;
10998  }
10999  if(*(u_begin) != value_type(0.0))
11000  {
11001  rho = slip::conj(*u_begin) / std::abs(*u_begin);
11002  }
11003  else
11004  {
11005  rho = value_type(1.0);
11006  }
11007  value_type scal = (rho/nu);
11008  slip::vector_scalar_multiplies(u_begin,u_end,scal,u_begin);
11009  *u_begin = value_type(1.0) + *u_begin;
11010  scal = value_type(1.0) / std::sqrt(std::real(*u_begin));
11011  slip::vector_scalar_multiplies(u_begin,u_end,scal,u_begin);
11012  *u_begin = std::real(*u_begin);
11013  nu = -(slip::conj(rho)) * nu;
11014 
11015  }
11016 
11017 
11047  template <typename RandomAccessIterator,
11048  typename MatrixIterator1>
11049  inline
11050  void
11051  right_householder_update(RandomAccessIterator V_first,
11052  RandomAccessIterator V_last,
11053  MatrixIterator1 M_up,
11054  MatrixIterator1 M_bot)
11055 
11056  {
11057  typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
11058 
11059  typename MatrixIterator1::difference_type sizeM = M_bot - M_up;
11060  std::size_t M_rows = sizeM[0];
11061 
11062 
11063  assert(sizeM[1] == (V_last - V_first));
11064 
11065  //computes w = MV
11066  slip::Array<value_type> w(M_rows);
11067  slip::matrix_vector_multiplies(M_up,M_bot,V_first,V_last,w.begin(),w.end());
11068 
11069  //computes M = M - beta w V^T
11070  slip::rank1_update(M_up,M_bot,
11071  value_type(-1.0),
11072  w.begin(),w.end(),
11073  V_first,V_last);
11074 
11075  }
11076 
11077 
11107  template <typename RandomAccessIterator,
11108  typename MatrixIterator1>
11109  inline
11110  void
11111  left_householder_update(RandomAccessIterator V_first,
11112  RandomAccessIterator V_last,
11113  MatrixIterator1 M_up,
11114  MatrixIterator1 M_bot)
11115 
11116  {
11117  typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
11118 
11119  typename MatrixIterator1::difference_type sizeM = M_bot - M_up;
11120 
11121  std::size_t M_cols = sizeM[1];
11122 
11123  assert(sizeM[0] == (V_last - V_first));
11124 
11125  //computes w = M^Hv
11126  slip::Array<value_type> w((std::size_t)sizeM[1]);
11127  for(std::size_t j = 0; j < M_cols; ++j)
11128  {
11129  w[j] = slip::hermitian_inner_product(M_up.col_begin(j),M_up.col_end(j),
11130  V_first,
11131  value_type(0));
11132  }
11133  //computes M = M - Vw^H
11134 
11135  slip::rank1_update(M_up,M_bot,
11136  value_type(-1.0),
11137  V_first,V_last,
11138  w.begin(),w.end());
11139  }
11140 
11151  template <typename Matrix1, typename Vector1,typename Matrix2>
11152  inline
11153  void
11154  left_householder_accumulate(const Matrix1& M,
11155  const Vector1& V0,
11156  Matrix2& Q)
11157 
11158  {
11159  typedef typename Matrix1::value_type value_type;
11160 
11161  std::size_t M_rows = M.rows();
11162  std::size_t M_cols = M.cols();
11163  int M_rows_m1 = int(M_rows) - 1;
11164  int M_cols_m1 = int(M_cols) - 1;
11165 
11166  slip::Array<value_type> V(M_rows);
11167  //init Q with the identity matrix
11168  slip::identity(Q);
11169  slip::Box2d<int> box;
11170  for(int j = M_rows-1; j > -1; --j)
11171  {
11172 
11173  V[j] = value_type(V0[j]);
11174  std::copy(M.col_begin(j)+(j+1),M.col_end(j),V.begin()+(j+1));
11175  box.set_coord(j,j,M_rows_m1,M_cols_m1);
11176  slip::left_householder_update(V.begin()+j,V.end(),
11177  Q.upper_left(box),Q.bottom_right(box));
11178 
11179  }
11180  }
11181 
11182 /* @} */
11183 
11185  /* @{ */
11209  template <typename Matrix1>
11210  inline
11211  void
11212  householder_hessenberg(Matrix1& M)
11213 
11214  {
11215  assert(M.rows() == M.cols());
11216  typedef typename Matrix1::value_type value_type;
11217 
11218  std::size_t M_rows = M.rows();
11219  std::size_t M_cols = M.cols();
11220 
11221  slip::Array<value_type> V(M_rows);
11222 
11223  slip::Box2d<int> box;
11224  slip::Box2d<int> box2;
11225  for(std::size_t j = 0; j < (M_cols-1); ++j)
11226  {
11227  //computes the householder vector from jth column of M
11228  box.set_coord(int(j+1),int(j),int(M_rows-1),int(M_cols-1));
11229  value_type beta = value_type(0);
11230 
11231  slip::housegen(M.upper_left(box).col_begin(0),
11232  M.upper_left(box).col_end(0),
11233  V.begin()+(j+1),V.end(),
11234  beta);
11235 
11236  //computes M = (I-VV^*)M(j+1:M_rows-1,j:M_cols-1)
11237  slip::left_householder_update(V.begin()+(j+1),V.end(),
11238  M.upper_left(box),M.bottom_right(box));
11239  box2.set_coord(int(0),int(j+1),int(M_rows-1),int(M_cols-1));
11240  //computes M = M(0:M_rows-1,j+1:M_cols-1)(I-betaVV^*)
11241  slip::right_householder_update(V.begin()+(j+1),V.end(),
11242  M.upper_left(box2),M.bottom_right(box2));
11243 
11244  }
11245  }
11246 
11247 
11276  template <typename Matrix1, typename Matrix2>
11277  inline
11278  void
11279  householder_hessenberg(const Matrix1& M, Matrix2& H)
11280 
11281  {
11282  assert(M.rows() == M.cols());
11283  assert(M.rows() == H.rows());
11284  assert(M.cols() == H.cols());
11285 
11286  typedef typename Matrix2::value_type value_type;
11287 
11288  //init H with M
11289  H = M;
11290  slip::householder_hessenberg(H);
11291 
11292  }
11293 
11294 
11334  template <typename Matrix1, typename Matrix2, typename Matrix3>
11335  inline
11336  void
11337  householder_hessenberg(const Matrix1& M, Matrix2& H, Matrix3& Q)
11338 
11339  {
11340  assert(M.rows() == M.cols());
11341  assert(H.rows() == H.cols());
11342  assert(Q.rows() == Q.cols());
11343  assert(M.rows() == H.rows());
11344  assert(M.rows() == Q.cols());
11345 
11346  typedef typename Matrix2::value_type value_type;
11347 
11348  //init Q with identity matrix
11349  slip::identity(Q);
11350 
11351  std::size_t M_rows = M.rows();
11352  std::size_t M_cols = M.cols();
11353 
11354  slip::Array<value_type> V(M_rows);
11355 
11356  slip::Box2d<int> box;
11357  slip::Box2d<int> box2;
11358  //init H with M
11359  std::copy(M.begin(),M.end(),H.begin());
11360  for(std::size_t j = 0; j < (M_cols-1); ++j)
11361  {
11362  //computes the householder vector from jth column of H
11363  box.set_coord(int(j+1),int(j),int(M_rows-1),int(M_cols-1));
11364  value_type beta = value_type(0);
11365 
11366  slip::housegen(H.upper_left(box).col_begin(0),
11367  H.upper_left(box).col_end(0),
11368  V.begin()+(j+1),
11369  V.end(),beta);
11370  //computes H = (I-VV^*)H(j+1:H_rows-1,j:H_cols-1)
11371  slip::left_householder_update(V.begin()+(j+1),V.end(),
11372  H.upper_left(box),H.bottom_right(box));
11373  std::fill(H.col_begin(j)+(j+2),H.col_end(j),value_type(0));
11374  box2.set_coord(int(0),int(j+1),int(M_rows-1),int(M_cols-1));
11375  //computes H = H(0:H_rows-1,j+1:H_cols-1)(I-VV^*)
11376  slip::right_householder_update(V.begin()+(j+1),V.end(),
11377  H.upper_left(box2),H.bottom_right(box2));
11378  //computes Q
11379  slip::right_householder_update(V.begin()+(j+1),V.end(),
11380  Q.upper_left(box2),Q.bottom_right(box2));
11381 
11382  }
11383  }
11384 
11385 /* @} */
11386 
11387 
11388 }//slip::
11389 
11390 
11391 #endif //SLIP_LINEAR_ALGEBRA_HPP
Array2d.hpp
Provides a class to manipulate 2d dynamic and generic arrays.
slip::hmatrix_matrix_multiplies
void hmatrix_matrix_multiplies(RandomAccessIterator2d1 M1_up, RandomAccessIterator2d1 M1_bot, RandomAccessIterator2d2 M2_up, RandomAccessIterator2d2 M2_bot, RandomAccessIterator2d3 result_up, RandomAccessIterator2d3 result_bot)
Computes the hermitian left multiplication of a matrix: .
Definition: linear_algebra.hpp:1289
slip::diagonal_solve
void diagonal_solve(RandomAccessIterator1 first, RandomAccessIterator1 last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision)
Solve the linear system D*X=B when D a diagonal matrix .
Definition: linear_algebra.hpp:6716
slip::is_skew_hermitian
bool is_skew_hermitian(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is skew hermitian.
Definition: linear_algebra.hpp:4084
slip::is_band_matrix
bool is_band_matrix(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename MatrixIterator1::size_type lower_band_width, const typename MatrixIterator1::size_type upper_band_width, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is a band matrix.
Definition: linear_algebra.hpp:4266
slip::is_complex
bool is_complex(const T &val)
Test if an element is complex.
Definition: linear_algebra_traits.hpp:179
slip::lower_band_solve
void lower_band_solve(MatrixIterator L_up, MatrixIterator L_bot, int p, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision)
Solve the linear system L*X=B when L is p lower banded .
Definition: linear_algebra.hpp:7273
slip::Box2d::set_coord
void set_coord(const CoordType &x11, const CoordType &x12, const CoordType &x21, const CoordType &x22)
Mutator of the coordinates of the Box2d.
Definition: Box2d.hpp:351
slip::cholesky_real
void cholesky_real(MatrixIterator1 A_up, MatrixIterator1 A_bot, MatrixIterator2 L_up, MatrixIterator2 L_bot, MatrixIterator3 LH_up, MatrixIterator3 LH_bot)
cholesky decomposition of a square real symmetric and positive definite Matrix. with L a lower trian...
Definition: linear_algebra.hpp:10008
slip::reduction
void reduction(Iterator1 seq1_beg, Iterator1 seq1_end, Iterator2 seq2_beg, U u)
Add u times the second sequence to the first one: seq1 = seq1 + u * seq2.
Definition: linear_algebra.hpp:2448
slip::DENSE_MATRIX_TYPE
DENSE_MATRIX_TYPE
Indicate the type of a dense matrix.
Definition: linear_algebra.hpp:134
slip::Array2d::upper_left
iterator2d upper_left()
Returns a read/write iterator2d that points to the first element of the Array2d. It points to the upp...
Definition: Array2d.hpp:2744
slip::imag
void imag(const Matrix1 &C, Matrix2 &I)
Extract the imaginary Matrix of a complex Matrix. .
Definition: linear_algebra.hpp:2684
slip::sign
T sign(T a)
Computes the sign of a.
Definition: macros.hpp:233
slip::Matrix::upper_left
iterator2d upper_left()
Returns a read/write iterator2d that points to the first element of the Matrix. It points to the uppe...
Definition: Matrix.hpp:3135
slip::right_householder_update
void right_householder_update(RandomAccessIterator V_first, RandomAccessIterator V_last, MatrixIterator1 M_up, MatrixIterator1 M_bot)
right multiplies the matrix M with the Householder matrix P:
Definition: linear_algebra.hpp:11051
slip::saxpy::operator()
_RT operator()(const _F &__x, const _S &__y) const
Definition: linear_algebra.hpp:171
slip::lu_solve
void lu_solve(const Matrix LU, const std::size_t nr, const std::size_t nc, const Vector1 &Indx, const std::size_t nv1, const Vector2 &B, const std::size_t nv2, Vector3 &X, const std::size_t nv3)
Solve Linear Equation using LU decomposition.
Definition: linear_algebra.hpp:5421
slip::lin_alg_traits::value_type
_Tp value_type
Definition: linear_algebra_traits.hpp:188
slip::band_lu
int band_lu(MatrixIterator1 A_up, MatrixIterator1 A_bot, const int p, const int q, RandomAccessIterator Ind_first, RandomAccessIterator Ind_last, typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type precision)
in place Band LU decomposition a square band Matrix A with L a p lower band matrix and U a q upper b...
Definition: linear_algebra.hpp:6773
slip::matrix_matrixt_multiplies
void matrix_matrixt_multiplies(RandomAccessIterator2d1 M1_up, RandomAccessIterator2d1 M1_bot, RandomAccessIterator2d2 Res_up, RandomAccessIterator2d2 Res_bot)
Computes the multiplication of a matrix with its transposate A^T .
Definition: linear_algebra.hpp:2047
slip::set_diagonal
void set_diagonal(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, Container2d &container, const int diag_number=0)
Set the diagonal [diag_first,diag_last) in the diagonal diag_number of a 2d container.
Definition: linear_algebra.hpp:3442
slip::Matrix::dim1
size_type dim1() const
Returns the number of rows (first dimension size) in the Matrix.
Definition: Matrix.hpp:3495
slip::lu_inv
void lu_inv(const Matrix1 &M, Matrix2 &IM)
Computes the inverse of a matrix using LU decomposition.
Definition: linear_algebra.hpp:5558
slip::iterator2d_box::x1
int x1() const
Access to the first index of the current iterator2d_box.
Definition: iterator2d_box.hpp:691
slip::Matrix::size_type
std::size_t size_type
Definition: Matrix.hpp:204
slip::givens_sin_cos
void givens_sin_cos(const Real &xi, const Real &xk, Real &sin, Real &cos)
Computes the Givens sinus and cosinus. .
Definition: linear_algebra.hpp:10720
slip::iterator2d_box::col_begin
Container2D::col_iterator col_begin(size_type col)
iterator2d_box element assignment operator.
Definition: iterator2d_box.hpp:579
slip::Vector::end
iterator end()
Returns a read/write iterator that points one past the last element in the Vector. Iteration is done in ordinary element order.
Definition: Vector.hpp:1481
slip::fill_diagonal
void fill_diagonal(Container2d &container, const typename Container2d::value_type &val, const int diag_number=0)
Fill the diagonal diag_number of a 2d container with the value val.
Definition: linear_algebra.hpp:3503
slip::L22_norm_cplx
slip::lin_alg_traits< typename std::iterator_traits< InputIterator >::value_type >::value_type L22_norm_cplx(InputIterator first, InputIterator last)
Computes the L22 norm of a complex container.
Definition: linear_algebra.hpp:3105
slip::Matrix::bottom_right
iterator2d bottom_right()
Returns a read/write iterator2d that points to the past the end element of the Matrix. It points to past the end element of the bottom right element of the Matrix.
Definition: Matrix.hpp:3151
slip::min
T & min(const GrayscaleImage< T > &M1)
Returns the min element of a GrayscaleImage.
Definition: GrayscaleImage.hpp:3693
linear_algebra_traits.hpp
slip::lower_bidiagonal_solve
void lower_bidiagonal_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
solve Ax = B with A lower bidiagonal
Definition: linear_algebra.hpp:8101
slip::is_null_diagonal
bool is_null_diagonal(MatrixIterator1 A_up, MatrixIterator1 A_bot, int diag_number, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix has a nul diagonal.
Definition: linear_algebra.hpp:4730
slip::LDLT_solve
void LDLT_solve(const Matrix1 &A, Vector1 &X, const Vector2 &B)
LDLT solve of system AX = B with A a square symmetric Matrix. with L a lower triangular matrix...
Definition: linear_algebra.hpp:9738
slip::unit_upper_bidiagonal_inv
void unit_upper_bidiagonal_inv(RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot)
Invert an unit upper bidiagonal matrix: .
Definition: linear_algebra.hpp:8690
slip::lu
int lu(const Matrix1 &M, Matrix2 &LU, Vector &Indx)
Computes the LU decomposition according to rows permutations of a matrix using Crout method LU is co...
Definition: linear_algebra.hpp:5164
slip::epsilon
slip::lin_alg_traits< T >::value_type epsilon()
Returns the epsilon value of a real or a complex.
Definition: linear_algebra_traits.hpp:211
slip::rotgen
void rotgen(Real &a, Real &b, Real &cos, Real &sin)
Computes the Givens sinus and cosinus.
Definition: linear_algebra.hpp:10664
slip::Array::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Array. Iteration is done in ord...
Definition: Array.hpp:1077
slip::cholesky_cplx
void cholesky_cplx(MatrixIterator1 A_up, MatrixIterator1 A_bot, MatrixIterator2 L_up, MatrixIterator2 L_bot, MatrixIterator3 LH_up, MatrixIterator3 LH_bot)
cholesky decomposition of a square hermitian positive definite Matrix. with L a lower triangular mat...
Definition: linear_algebra.hpp:9822
slip::Array2d::bottom_right
iterator2d bottom_right()
Returns a read/write iterator2d that points to the past the end element of the Array2d. It points to past the end element of the bottom right element of the Array2d.
Definition: Array2d.hpp:2759
slip::is_lower_bidiagonal
bool is_lower_bidiagonal(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is lower_bidiagonal.
Definition: linear_algebra.hpp:4470
slip::is_squared
bool is_squared(const Matrix &A)
Test if a matrix is squared.
Definition: linear_algebra.hpp:3800
slip::Matrix::cols
size_type cols() const
Returns the number of columns (second dimension size) in the Matrix.
Definition: Matrix.hpp:3512
slip::is_tridiagonal
bool is_tridiagonal(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is tridiagonal.
Definition: linear_algebra.hpp:4364
slip::analyse_matrix
slip::DENSE_MATRIX_TYPE analyse_matrix(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Analyse the type of a matrix.
Definition: linear_algebra.hpp:4785
slip::complex_sign
Definition: macros.hpp:258
slip::conjz1z2
Computes the complex operation: ( ) between two complex values z1 and z2.
Definition: linear_algebra.hpp:218
slip::z1conjz2::operator()
_Tp operator()(const _Tp &z1, const _Tp &z2) const
Definition: linear_algebra.hpp:198
slip::is_upper_bidiagonal
bool is_upper_bidiagonal(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is upper_bidiagonal.
Definition: linear_algebra.hpp:4418
slip::hmatrix_vector_multiplies
void hmatrix_vector_multiplies(RandomAccessIterator2d M_up, RandomAccessIterator2d M_bot, RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, RandomAccessIterator2 R_first, RandomAccessIterator2 R_last)
Computes the hermitian matrix left multiplication of a vector: .
Definition: linear_algebra.hpp:1403
slip::is_upper_triangular
bool is_upper_triangular(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is upper_triangular.
Definition: linear_algebra.hpp:4628
slip::matrix_vector_multiplies
void matrix_vector_multiplies(RandomAccessIterator2d M_upper_left, RandomAccessIterator2d M_bottom_right, RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 result_first, RandomAccessIterator2 result_last)
Computes the multiplication of a Matrix and a Vector: Result=MV.
Definition: linear_algebra.hpp:963
slip::upper_triangular_inv
void upper_triangular_inv(MatrixIterator1 A_up, MatrixIterator1 A_bot, MatrixIterator2 Ainv_up, MatrixIterator2 Ainv_bot, typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type precision=typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type(1.0E-6))
Invert an upper triangular Matrix .
Definition: linear_algebra.hpp:5895
slip::lower_triangular_inv
void lower_triangular_inv(MatrixIterator1 A_up, MatrixIterator1 A_bot, MatrixIterator2 Ainv_up, MatrixIterator2 Ainv_bot, typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type precision=typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type(1.0E-6))
Invert a lower triangular Matrix .
Definition: linear_algebra.hpp:6212
slip::thomas_inv
void thomas_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
Invert a tridiagonal squared matrix A with the Thomas method .
Definition: linear_algebra.hpp:9000
slip::upper_triangular_solve
void upper_triangular_solve(MatrixIterator U_up, MatrixIterator U_bot, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision)
Solve the linear system U*X=B when U is upper triangular .
Definition: linear_algebra.hpp:5722
slip::iota
void iota(ForwardIterator first, ForwardIterator last, T value, T step=T(1))
Iota assigns sequential increasing values based on a predefined step to a range. That is...
Definition: arithmetic_op.hpp:2087
slip::thomas_solve
void thomas_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
Solve the tridiagonal system A*X=B with the Thomas method .
Definition: linear_algebra.hpp:8795
slip::Matrix::col_end
col_iterator col_end(const size_type col)
Returns a read/write iterator that points one past the end element of the column column in the Matrix...
slip::is_upper_hessenberg
bool is_upper_hessenberg(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is upper_hessenberg.
Definition: linear_algebra.hpp:4521
slip::conjz1z2::operator()
_Tp operator()(const _Tp &z1, const _Tp &z2) const
Definition: linear_algebra.hpp:222
slip::real
void real(const Matrix1 &C, Matrix2 &R)
Extract the real Matrix of a complex Matrix. .
Definition: linear_algebra.hpp:2643
slip::abs
HyperVolume< T > abs(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5064
slip::lower_triangular_solve
void lower_triangular_solve(MatrixIterator L_up, MatrixIterator L_bot, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision)
Solve the linear system L*X=B when L is lower triangular .
Definition: linear_algebra.hpp:6048
slip::matrix_shift
void matrix_shift(RandomAccessIterator2d1 M_upper_left, RandomAccessIterator2d1 M_bottom_right, const T &lambda, RandomAccessIterator2d2 R_upper_left, RandomAccessIterator2d2 R_bottom_right)
Computes with In the identity matrix.
Definition: linear_algebra.hpp:2146
slip::sin
HyperVolume< T > sin(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5101
slip::unit_upper_triangular_solve
void unit_upper_triangular_solve(MatrixIterator U_up, MatrixIterator U_bot, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last)
Solve the linear system U*X=B when U is unit upper triangular .
Definition: linear_algebra.hpp:6363
slip::left_householder_update
void left_householder_update(RandomAccessIterator V_first, RandomAccessIterator V_last, MatrixIterator1 M_up, MatrixIterator1 M_bot)
Left multiplies the Householder matrix P with the matrix M: .
Definition: linear_algebra.hpp:11111
Array.hpp
Provides a class to manipulate 1d dynamic and generic arrays.
slip::cos
HyperVolume< T > cos(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5083
slip::partial_pivot
int partial_pivot(MatrixIterator M_up, MatrixIterator M_bot)
Returns the partial pivot of a 2d range.
Definition: linear_algebra.hpp:4909
slip::unit_upper_band_solve
void unit_upper_band_solve(MatrixIterator U_up, MatrixIterator U_bot, int q, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision)
Solve the linear system U*X=B when U is unit q width upper banded .
Definition: linear_algebra.hpp:7154
slip::complex_left_givens
void complex_left_givens(Matrix &M, const SizeType &row1, const SizeType &row2, const Complex &sinus, const typename Complex::value_type &cosinus, const SizeType &col1, const SizeType &col2)
Apply complex left Givens rotation multiplication.
Definition: linear_algebra.hpp:10885
slip::rank1_update
void rank1_update(RandomAccessIterator2d1 A_up, RandomAccessIterator2d1 A_bot, const T &alpha, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 Y_first, RandomAccessIterator2 Y_last)
Computes .
Definition: linear_algebra.hpp:2346
slip::saxpy::saxpy
saxpy(const _AT &a)
Definition: linear_algebra.hpp:166
slip::tridiagonal_cholesky_solve
void tridiagonal_cholesky_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last)
Solve the tridiagonal symmetric positive definite system T*X=B with the Cholesky method where: is t...
Definition: linear_algebra.hpp:9321
slip::gen_matrix_vector_multiplies
void gen_matrix_vector_multiplies(RandomAccessIterator2d M_up, RandomAccessIterator2d M_bot, RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, RandomAccessIterator2 Y_first, RandomAccessIterator2 Y_last)
Computes the generalised multiplication of a Matrix and a Vector: computes the update ...
Definition: linear_algebra.hpp:1510
slip::is_diagonal
bool is_diagonal(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is diagonal.
Definition: linear_algebra.hpp:3899
slip::SINGULAR_MATRIX_ERROR
const std::string SINGULAR_MATRIX_ERROR
Definition: error.hpp:84
slip::right_givens
void right_givens(Matrix &M, const SizeType &col1, const SizeType &col2, const Real &sinus, const Real &cosinus, const SizeType &row1, const SizeType &row2)
Apply right Givens rotation multiplication.
Definition: linear_algebra.hpp:10841
slip::left_givens
void left_givens(Matrix &M, const SizeType &row1, const SizeType &row2, const Real &sinus, const Real &cosinus, const SizeType &col1, const SizeType &col2)
Apply left Givens rotation multiplication.
Definition: linear_algebra.hpp:10798
slip::Array::end
iterator end()
Returns a read/write iterator that points one past the last element in the Array. Iteration is done i...
Definition: Array.hpp:1084
slip::upper_bidiagonal_inv
void upper_bidiagonal_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
Invert an upper bidiagonal squared matrix .
Definition: linear_algebra.hpp:8451
slip::gauss_solve_with_filling
bool gauss_solve_with_filling(const Matrix &M, Vector1 &X, const Vector2 &B, typename slip::lin_alg_traits< typename Matrix::value_type >::value_type precision)
Solve the linear system M*X=B with the Gauss pivot method, the null pivot are replaced by precision...
Definition: linear_algebra.hpp:7687
slip::Vector::size
size_type size() const
Returns the number of elements in the Vector.
Definition: Vector.hpp:1743
slip::row_norm
slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator2d >::value_type >::value_type row_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
Computes the row norm ( ) of a 2d range.
Definition: linear_algebra.hpp:3141
slip::upper_bidiagonal_solve
void upper_bidiagonal_solve(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
solve Ax = B with A upper bidiagonal
Definition: linear_algebra.hpp:8360
slip::multiplies
void multiplies(InputIterator1 __first1, InputIterator1 __last1, InputIterator2 __first2, OutputIterator __result)
Computes the pointwise product of two ranges.
Definition: arithmetic_op.hpp:508
slip::lower_bidiagonal_inv
void lower_bidiagonal_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
Invert a lower bidiagonal matrix: .
Definition: linear_algebra.hpp:8182
slip::vector_scalar_multiplies
void vector_scalar_multiplies(RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, const T &scal, RandomAccessIterator2 result_first)
Computes the multiplication of a vector V by a scalar scal.
Definition: linear_algebra.hpp:688
slip::pivot_max
int pivot_max(MatrixIterator M_up, MatrixIterator M_bot)
Returns the row position of the maximum partial pivot of a 2d range.
Definition: linear_algebra.hpp:4877
slip::copy
void copy(_II first, _II last, _OI output_first)
Copy algorithm optimized for slip iterators.
Definition: copy_ext.hpp:177
slip::DPoint2d
Difference of Point2D class, specialization of DPoint<CoordType,DIM> with DIM = 2.
Definition: Array2d.hpp:129
complex_cast.hpp
Provides some macros which are used for using complex as real.
slip::cholesky_inv
void cholesky_inv(MatrixIterator1 A_up, MatrixIterator1 A_bot, MatrixIterator2 Ainv_up, MatrixIterator2 Ainv_bot, typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type precision=typename slip::lin_alg_traits< typename MatrixIterator1::value_type >::value_type(1.0E-6))
cholesky inverse of a square hermitian positive definite Matrix A.
Definition: linear_algebra.hpp:10550
slip::Vector
Numerical Vector class. This container statisfies the RandomAccessContainer concepts of the Standard ...
Definition: Vector.hpp:104
slip::frobenius_norm
slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator2d >::value_type >::value_type frobenius_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
Computes the Frobenius norm ( ) of a 2d range.
Definition: linear_algebra.hpp:3286
slip::z1conjz2
Computes the complex operation: ( ) between two complex values z1 and z2.
Definition: linear_algebra.hpp:194
slip::cholesky_solve
void cholesky_solve(MatrixIterator A_up, MatrixIterator A_bot, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision=typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type(1.0E-6))
cholesky solve of system AX = B with A a square hermitian positive definite Matrix.
Definition: linear_algebra.hpp:10412
slip::is_lower_triangular
bool is_lower_triangular(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is lower_triangular.
Definition: linear_algebra.hpp:4679
slip::DPoint2d::set_coord
void set_coord(const CoordType &dx1, const CoordType &dx2)
Accessor/Mutator of the coordinates of DPoint2d.
Definition: DPoint2d.hpp:253
slip::housegen
void housegen(VectorIterator1 a_begin, VectorIterator1 a_end, VectorIterator2 u_begin, VectorIterator2 u_end, typename std::iterator_traits< VectorIterator1 >::value_type &nu)
Compute the Householder vector u of a vector a.
Definition: linear_algebra.hpp:10981
macros.hpp
Provides some mathematical functors and constants.
slip::sqrt
HyperVolume< T > sqrt(const HyperVolume< T > &M)
Definition: HyperVolume.hpp:5074
slip::unit_lower_triangular_solve
void unit_lower_triangular_solve(MatrixIterator L_up, MatrixIterator L_bot, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last)
Solve the linear system L*X=B when L is unit lower triangular .
Definition: linear_algebra.hpp:6539
slip::identity
Matrix identity(const std::size_t nr, const std::size_t nc)
Returns an identity matrix which dimensions are nr*nc.
Definition: linear_algebra.hpp:3608
slip::unit_lower_bidiagonal_solve
void unit_lower_bidiagonal_solve(RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last)
solve Ax = B with A unit lower bidiagonal
Definition: linear_algebra.hpp:7946
slip::matrix_scalar_multiplies
void matrix_scalar_multiplies(RandomAccessIterator2d1 M_upper_left, RandomAccessIterator2d1 M_bottom_right, const T &scal, RandomAccessIterator2d2 R_upper_left, RandomAccessIterator2d2 R_bottom_right)
Computes the multiplication of a Matrix by a scalar R = M*scal.
Definition: linear_algebra.hpp:856
slip::gen_matrix_matrix_multiplies
void gen_matrix_matrix_multiplies(RandomAccessIterator2d1 A_up, RandomAccessIterator2d1 A_bot, RandomAccessIterator2d2 B_up, RandomAccessIterator2d2 B_bot, RandomAccessIterator2d3 C_up, RandomAccessIterator2d3 C_bot)
Computes the generalized multiplication of a two Matrix: .
Definition: linear_algebra.hpp:1741
slip::Matrix::row_end
row_iterator row_end(const size_type row)
Returns a read/write iterator that points one past the end element of the row row in the Matrix...
slip::matrix_matrix_multiplies
void matrix_matrix_multiplies(MatrixIterator1 M1_up, MatrixIterator1 M1_bot, MatrixIterator2 M2_up, MatrixIterator2 M2_bot, MatrixIterator3 R_up, MatrixIterator3 R_bot)
Computes the matrix matrix multiplication of 2 2d ranges.
Definition: linear_algebra.hpp:1621
slip::is_identity
bool is_identity(const Matrix &A, const typename slip::lin_alg_traits< typename Matrix::value_type >::value_type tol=slip::epsilon< typename Matrix::value_type >())
Test if a matrix is identity.
Definition: linear_algebra.hpp:3848
slip::Matrix
Numerical matrix class. This container statisfies the BidirectionnalContainer concepts of the STL...
Definition: GenericMultiComponent2d.hpp:132
slip::unit_upper_bidiagonal_solve
void unit_upper_bidiagonal_solve(RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2 X_first, RandomAccessIterator2 X_last, RandomAccessIterator3 B_first, RandomAccessIterator3 B_last)
solve Ax = B with A unit upper bidiagonal
Definition: linear_algebra.hpp:8620
arithmetic_op.hpp
Provides some algorithms to computes arithmetical operations on ranges.
error.hpp
Provides SLIP error messages.
slip::gauss_solve
bool gauss_solve(const Matrix &M, Vector1 &X, const Vector2 &B, typename slip::lin_alg_traits< typename Matrix::value_type >::value_type precision)
Solve the linear system M*X=B using the Gauss elemination partial pivoting.
Definition: linear_algebra.hpp:7566
slip::left_householder_accumulate
void left_householder_accumulate(const Matrix1 &M, const Vector1 &V0, Matrix2 &Q)
Computes Q = Q1Q2...Qn from the inplace Householder QR decomposition.
Definition: linear_algebra.hpp:11154
slip::tridiagonal_cholesky_inv
void tridiagonal_cholesky_inv(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot)
Invert the tridiagonal symmetric positive definite system T: where: is the conjugate of and ...
Definition: linear_algebra.hpp:9516
slip::POSITIVE_DEFINITE_MATRIX_ERROR
const std::string POSITIVE_DEFINITE_MATRIX_ERROR
Definition: error.hpp:85
slip::lu_det
Matrix1::value_type lu_det(const Matrix1 &M)
Computes the determinant of a matrix using LU decomposition.
Definition: linear_algebra.hpp:5619
slip::inverse
void inverse(const Matrix1 &M, const std::size_t nr1, const std::size_t nc1, Matrix2 &IM, const std::size_t nr2, const std::size_t nc2)
Computes the inverse of a matrix using gaussian elimination.
Definition: linear_algebra.hpp:4963
slip::Matrix::rows
size_type rows() const
Returns the number of rows (first dimension size) in the Matrix.
Definition: Matrix.hpp:3499
slip::hermitian_inner_product
std::iterator_traits< RandomAccessIterator1 >::value_type hermitian_inner_product(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 first2, typename std::iterator_traits< RandomAccessIterator1 >::value_type init)
Computes the hermitian inner-product of two ranges X and Y: .
Definition: linear_algebra.hpp:271
slip::tridiagonal_lu
void tridiagonal_lu(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2 L_low_first, RandomAccessIterator2 L_low_last, RandomAccessIterator3 U_up_first, RandomAccessIterator3 U_up_last, typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type precision=typename slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator1 >::value_type >::value_type(1.0E-6))
Computes the LU decomposition for a tridiagonal matrix .
Definition: linear_algebra.hpp:7846
slip::Matrix::col_begin
col_iterator col_begin(const size_type col)
Returns a read/write iterator that points to the first element of the column column in the Matrix...
slip::upper_band_solve
void upper_band_solve(MatrixIterator U_up, MatrixIterator U_bot, int q, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, typename slip::lin_alg_traits< typename MatrixIterator::value_type >::value_type precision)
Solve the linear system U*X=B when U is q width upper banded .
Definition: linear_algebra.hpp:7025
slip::householder_hessenberg
void householder_hessenberg(Matrix1 &M)
Householder Hessenberg reduction of the square matrix M. The result is overwritten in M...
Definition: linear_algebra.hpp:11212
slip::inner_product
std::iterator_traits< RandomAccessIterator1 >::value_type inner_product(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 first2, typename std::iterator_traits< RandomAccessIterator1 >::value_type init=typename std::iterator_traits< RandomAccessIterator1 >::value_type())
Computes the inner_product of two ranges X and Y: .
Definition: linear_algebra.hpp:323
slip::get_diagonal
void get_diagonal(const Container2d &container, RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, const int diag_number=0)
Get the diagonal diag_number of a 2d container.
Definition: linear_algebra.hpp:3370
slip::hvector_matrix_multiplies
void hvector_matrix_multiplies(RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, RandomAccessIterator2d M_up, RandomAccessIterator2d M_bot, RandomAccessIterator2 result_first, RandomAccessIterator2 result_last)
Computes the hermitian left multiplication: .
Definition: linear_algebra.hpp:1078
slip::is_symmetric
bool is_symmetric(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is symmetric.
Definition: linear_algebra.hpp:4189
slip::is_lower_hessenberg
bool is_lower_hessenberg(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is lower_hessenberg.
Definition: linear_algebra.hpp:4575
slip::gen_inner_product
std::iterator_traits< RandomAccessIterator2 >::value_type gen_inner_product(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2d A_upper_left, RandomAccessIterator2d A_bottom_right, RandomAccessIterator2 first2, RandomAccessIterator2 last2, typename std::iterator_traits< RandomAccessIterator2 >::value_type init=typename std::iterator_traits< RandomAccessIterator2 >::value_type())
Computes the generalized inner_product of two ranges: .
Definition: linear_algebra.hpp:442
slip::tvector_matrix_multiplies
void tvector_matrix_multiplies(RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, RandomAccessIterator2d M_up, RandomAccessIterator2d M_bot, RandomAccessIterator2 result_first, RandomAccessIterator2 result_last)
Computes the transpose left multiplication: .
Definition: linear_algebra.hpp:1189
slip::rayleigh
std::iterator_traits< RandomAccessIterator1 >::value_type rayleigh(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2d A_upper_left, RandomAccessIterator2d A_bottom_right)
Computes the rayleigh coefficient of a vector x: .
Definition: linear_algebra.hpp:557
slip::inner_product
Vector1::value_type inner_product(const Vector1 &V1, const Vector2 &V2)
Computes the inner_product of two ranges X and Y: .
Definition: linear_algebra.hpp:381
slip::cholesky
void cholesky(const Matrix1 &A, Matrix2 &L, Matrix3 &LT)
cholesky decomposition of a square hermitian positive definite Matrix. with L a lower triangular mat...
Definition: linear_algebra.hpp:10290
apply.hpp
Provides some algorithms to apply C like functions to ranges.
slip::conj_vector_scalar_multiplies
void conj_vector_scalar_multiplies(RandomAccessIterator1 V_first, RandomAccessIterator1 V_last, const T &scal, RandomAccessIterator2 result_first)
Computes the multiplication of a vector by a scalar scal.
Definition: linear_algebra.hpp:774
slip::transpose
void transpose(const Matrix1 &M, const std::size_t nr1, const std::size_t nc1, Matrix2 &TM, const std::size_t nr2, const std::size_t nc2)
Computes the transpose of a matrix .
Definition: linear_algebra.hpp:2516
slip::LDLT_decomposition
void LDLT_decomposition(Container2d &A)
in place LU decomposition for symmetric matrix
Definition: linear_algebra.hpp:9633
slip::col_norm
slip::lin_alg_traits< typename std::iterator_traits< RandomAccessIterator2d >::value_type >::value_type col_norm(RandomAccessIterator2d upper_left, RandomAccessIterator2d bottom_right)
Computes the column norm ( ) of a 2d range.
Definition: linear_algebra.hpp:3214
slip::complex_right_givens
void complex_right_givens(Matrix &M, const SizeType &col1, const SizeType &col2, const Complex &sinus, const typename Complex::value_type &cosinus, const SizeType &row1, const SizeType &row2)
Apply complex right Givens rotation multiplication.
Definition: linear_algebra.hpp:10928
slip::rank1_tensorial_product
void rank1_tensorial_product(RandomAccessIterator1 base1_first, RandomAccessIterator1 base1_end, RandomAccessIterator2 base2_first, RandomAccessIterator2 base2_end, RandomAccessIterator2d matrix_upper_left, RandomAccessIterator2d matrix_bottom_right)
Computes the tensorial product of two rank one tensors it provides a rank2 tensor (Container2d) of s...
Definition: linear_algebra.hpp:2859
slip::Array2d::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Array2d. Iteration is done in o...
Definition: Array2d.hpp:2345
slip::Matrix::dim2
size_type dim2() const
Returns the number of columns (second dimension size) in the Matrix.
Definition: Matrix.hpp:3504
slip::tridiagonal_cholesky
void tridiagonal_cholesky(RandomAccessIterator1 diag_first, RandomAccessIterator1 diag_last, RandomAccessIterator1 up_diag_first, RandomAccessIterator1 up_diag_last, RandomAccessIterator2 R_diag_first, RandomAccessIterator2 R_diag_last, RandomAccessIterator3 R_low_first, RandomAccessIterator3 R_low_last)
Tridiagonal symmetric positive decomposition where: is the conjugate of and .
Definition: linear_algebra.hpp:9203
slip::hermitian_transpose
void hermitian_transpose(const Matrix1 &M, Matrix2 &TM)
Computes the hermitian transpose of a matrix which is the transpose of the conjugate.
Definition: linear_algebra.hpp:2476
slip::unit_lower_band_solve
void unit_lower_band_solve(MatrixIterator L_up, MatrixIterator L_bot, int p, RandomAccessIterator1 X_first, RandomAccessIterator1 X_last, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last)
Solve the linear system L*X=B when L is unit p lower banded .
Definition: linear_algebra.hpp:7400
slip::Array
This is a linear (one-dimensional) dynamic template container. This container statisfies the RandomAc...
Definition: Array.hpp:94
slip::is_hermitian
bool is_hermitian(MatrixIterator1 A_up, MatrixIterator1 A_bot, const typename slip::lin_alg_traits< typename std::iterator_traits< MatrixIterator1 >::value_type >::value_type tol=slip::epsilon< typename std::iterator_traits< MatrixIterator1 >::value_type >())
Test if a matrix is hermitian.
Definition: linear_algebra.hpp:3986
slip::hilbert
void hilbert(Matrix &A)
Replaces A with the Hilbert matrix: .
Definition: linear_algebra.hpp:3717
slip::aXpY
void aXpY(const AT &a, RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 first2)
Computes the saxpy ("scalar a x plus b") for each element of two ranges x and y.
Definition: linear_algebra.hpp:632
slip::Matrix::row_begin
row_iterator row_begin(const size_type row)
Returns a read/write iterator that points to the first element of the row row in the Matrix...
slip::toeplitz
void toeplitz(RandomAccessIterator first, RandomAccessIterator last, Matrix &A)
Constructs the Toeplitz matrix from given a range.
Definition: linear_algebra.hpp:3758
slip::complex_givens_sin_cos
void complex_givens_sin_cos(const Complex &xi, const Complex &xk, Complex &sin, typename Complex::value_type &cos)
Computes the complex Givens sinus and cosinus. with .
Definition: linear_algebra.hpp:10761
slip::lin_alg_traits
Definition: linear_algebra_traits.hpp:186
slip::conj
int conj(const int arg)
Definition: complex_cast.hpp:122
slip::apply
void apply(InputIterator first, InputIterator last, OutputIterator result, typename std::iterator_traits< OutputIterator >::value_type(*function)(typename std::iterator_traits< InputIterator >::value_type))
Applies a C-function to each element of a range.
Definition: apply.hpp:128
slip::Array2d
This is a two-dimensional dynamic and generic container. This container statisfies the Bidirectionnal...
Definition: Array2d.hpp:135
slip::max
T & max(const GrayscaleImage< T > &M1)
Returns the max element of a GrayscaleImage.
Definition: GrayscaleImage.hpp:3704
slip::unit_lower_bidiagonal_inv
void unit_lower_bidiagonal_inv(RandomAccessIterator1 low_diag_first, RandomAccessIterator1 low_diag_last, RandomAccessIterator2d Ainv_up, RandomAccessIterator2d Ainv_bot)
Invert a lower unit bidiagonal matrix: .
Definition: linear_algebra.hpp:8007
slip::Vector::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Vector. Iteration is done in or...
Definition: Vector.hpp:1474
slip::outer_product
void outer_product(VectorIterator1 V_begin, VectorIterator1 V_end, VectorIterator2 W_begin, VectorIterator2 W_end, MatrixIterator R_up, MatrixIterator R_bot)
Computes the hermitian outer product of the vector V and the transpose of the conjugate of the vector...
Definition: linear_algebra.hpp:2739
slip::saxpy
Computes the saxpy ("scalar a x plus b") operation: ( ) between two values.
Definition: linear_algebra.hpp:162
norms.hpp
Provides some algorithms to compute norms.
Generated on Tue Sep 30 2014 10:05:57 for SLIP by   doxygen 1.8.6