23 #ifndef O2SCL_VECTOR_H 24 #define O2SCL_VECTOR_H 57 #include <gsl/gsl_vector.h> 58 #include <gsl/gsl_sys.h> 59 #include <gsl/gsl_matrix.h> 61 #include <o2scl/misc.h> 62 #include <o2scl/uniform_grid.h> 76 d=(
const double *)m->data;
79 double operator[](
size_t i)
const {
99 d=(
const double *)m->data;
104 double operator()(
size_t i,
size_t j)
const {
126 template<
class vec_t,
class vec2_t>
129 if (dest.size()<N) dest.resize(N);
134 for(i=m;i+3<N;i+=4) {
155 template<
class vec_t,
class vec2_t>
161 for(i=m;i+3<N;i+=4) {
178 template<
class mat_t,
class mat2_t>
180 size_t m=src.size1();
181 size_t n=src.size2();
182 if (dest.size1()<m || dest.size2()<n) dest.resize(m,n);
183 for(
size_t i=0;i<m;i++) {
184 for(
size_t j=0;j<n;j++) {
203 template<
class mat_t,
class mat2_t>
205 for(
size_t i=0;i<M;i++) {
206 for(
size_t j=0;j<N;j++) {
224 template<
class mat_t,
class mat2_t>
226 size_t m=src.size1();
227 size_t n=src.size2();
228 if (dest.size1()<n || dest.size2()<m) dest.resize(n,m);
229 for(
size_t i=0;i<m;i++) {
230 for(
size_t j=0;j<n;j++) {
245 template<
class mat_t,
class mat2_t>
247 for(
size_t i=0;i<m;i++) {
248 for(
size_t j=0;j<n;j++) {
264 template<
class mat_t,
class data_t>
266 size_t m=src.size1();
267 size_t n=src.size2();
271 for(
size_t i=0;i<n;i++) {
272 for(
size_t j=0;j<n;j++) {
289 template<
class mat_t,
class data_t>
292 for(
size_t i=0;i<m;i++) {
293 for(
size_t j=0;j<n;j++) {
307 size_t m=src.size1();
308 size_t n=src.size2();
310 for(
size_t i=0;i<m;i++) {
311 for(
size_t j=i+1;j<n;j++) {
312 if (src(i,j)!=0.0)
return false;
321 size_t m=src.size1();
322 size_t n=src.size2();
324 for(
size_t j=0;j<n;j++) {
325 for(
size_t i=j+1;i<m;i++) {
326 if (src(i,j)!=0.0)
return false;
336 size_t m=src.size1();
337 size_t n=src.size2();
338 for(
size_t i=0;i<m;i++) {
339 for(
size_t j=i+1;j<n;j++) {
350 size_t m=src.size1();
351 size_t n=src.size2();
352 for(
size_t j=0;j<n;j++) {
353 for(
size_t i=j+1;i<m;i++) {
366 for(
size_t i=0;i<m;i++) {
367 for(
size_t j=i+1;j<n;j++) {
368 if (src(i,j)!=0.0)
return false;
380 for(
size_t j=0;j<n;j++) {
381 for(
size_t i=j+1;i<m;i++) {
382 if (src(i,j)!=0.0)
return false;
393 for(
size_t i=0;i<m;i++) {
394 for(
size_t j=i+1;j<n;j++) {
406 for(
size_t j=0;j<n;j++) {
407 for(
size_t i=j+1;i<m;i++) {
422 template<
class vec_t,
class vec2_t,
class data_t>
431 for(i=m;i+3<N;i+=4) {
458 template<
class vec_t,
class vec2_t,
class data_t>
461 if (v2.size()<N) N=v2.size();
469 for(i=m;i+3<N;i+=4) {
492 template<
class vec_t,
class vec2_t>
494 return vector_swap<vec_t,vec2_t,double>(N,v1,v2);
508 template<
class vec_t,
class vec2_t>
510 return vector_swap<vec_t,vec2_t,double>(v1,v2);
518 template<
class vec_t,
class data_t>
533 template<
class vec_t>
535 return vector_swap<vec_t,double>(v,i,j);
544 template<
class mat_t,
class mat2_t,
class data_t>
547 for(
size_t i=0;i<M;i++) {
548 for(
size_t j=0;j<N;j++) {
563 template<
class mat_t,
class mat2_t,
class data_t>
565 return matrix_swap<mat_t,mat2_t,double>(M,N,m1,m2);
573 template<
class mat_t,
class data_t>
574 void matrix_swap(mat_t &m,
size_t i1,
size_t j1,
size_t i2,
size_t j2) {
575 data_t temp=m(i1,j1);
586 template<
class mat_t>
588 size_t i2,
size_t j2) {
589 return matrix_swap<mat_t,double>(m,i1,j1,i2,j2);
597 template<
class mat_t,
class data_t>
600 for(
size_t i=0;i<M;i++) {
614 template<
class mat_t>
616 return matrix_swap_cols<mat_t,double>(M,m,j1,j2);
625 template<
class mat_t,
class data_t>
628 for(
size_t j=0;j<N;j++) {
642 template<
class mat_t>
644 return matrix_swap_rows<mat_t,double>(N,m,i1,i2);
652 template<
class vec_t,
class data_t>
660 if (j<n && data[j] < data[j+1]) j++;
661 if (!(v < data[j]))
break;
707 template<
class vec_t,
class data_t>
720 sort_downheap<vec_t,data_t>(data,N,k);
728 sort_downheap<vec_t,data_t>(data,N,0);
736 template<
class vec_t,
class vec_
size_t>
740 const size_t pki = order[k];
745 if (j < N && data[order[j]] < data[order[j + 1]]) {
750 if (!(data[pki] < data[order[j]])) {
802 template<
class vec_t,
class vec_
size_t>
811 for (i = 0 ; i < n ; i++) {
825 sort_index_downheap<vec_t,vec_size_t>(N,data,order,k);
831 size_t tmp = order[0];
838 sort_index_downheap<vec_t,vec_size_t>(N,data,order,0);
882 template<
class vec_t,
class vec_
size_t>
897 template<
class vec_t>
899 return vector_sort<vec_t,double>(n,data);
922 template<
class vec_t,
class data_t>
925 O2SCL_ERR2(
"Subset length greater than size in ",
935 data_t xbound=data[0];
939 for(
size_t i=1;i<n;i++) {
943 }
else if (xi>=xbound) {
947 for(i1=j-1;i1>0;i1--) {
948 if (xi>smallest[i1-1])
break;
949 smallest[i1]=smallest[i1-1];
952 xbound=smallest[j-1];
973 template<
class vec_t,
class data_t>
975 size_t n=data.size();
976 if (smallest.size()<k) smallest.resize(k);
995 template<
class vec_t,
class data_t,
class vec_
size_t>
999 O2SCL_ERR2(
"Subset length greater than size in ",
1000 "function vector_smallest_index().",
exc_einval);
1004 "function vector_smallest_index().",
exc_einval);
1011 data_t xbound=data[0];
1015 for(
size_t i=1;i<n;i++) {
1019 }
else if (xi>=xbound) {
1023 for(i1=j-1;i1>0;i1--) {
1024 if (xi>data[index[i1-1]])
break;
1025 index[i1]=index[i1-1];
1028 xbound=data[index[j-1]];
1035 template<
class vec_t,
class data_t,
class vec_
size_t>
1037 vec_size_t &index) {
1038 size_t n=data.size();
1039 if (index.size()<k) index.resize(k);
1040 return o2scl::vector_smallest_index<vec_t,data_t,vec_size_t>
1060 template<
class vec_t,
class data_t>
1063 O2SCL_ERR2(
"Subset length greater than size in ",
1073 data_t xbound=data[0];
1077 for(
size_t i=1;i<n;i++) {
1081 }
else if (xi<=xbound) {
1085 for(i1=j-1;i1>0;i1--) {
1086 if (xi<largest[i1-1])
break;
1087 largest[i1]=largest[i1-1];
1090 xbound=largest[j-1];
1111 template<
class vec_t,
class data_t>
1113 size_t n=data.size();
1114 if (largest.size()<k) largest.resize(k);
1123 template<
class vec_t,
class data_t>
1130 for(
size_t i=1;i<n;i++) {
1140 template<
class vec_t,
class data_t>
1143 size_t n=data.size();
1148 for(
size_t i=1;i<n;i++) {
1159 template<
class vec_t,
class data_t>
1167 for(
size_t i=1;i<n;i++) {
1178 template<
class vec_t,
class data_t>
1187 for(
size_t i=1;i<n;i++) {
1198 template<
class vec_t,
class data_t>
1205 for(
size_t i=1;i<n;i++) {
1215 template<
class vec_t,
class data_t>
1218 size_t n=data.size();
1223 for(
size_t i=1;i<n;i++) {
1234 template<
class vec_t,
class data_t>
1242 for(
size_t i=1;i<n;i++) {
1253 template<
class vec_t,
class data_t>
1262 for(
size_t i=1;i<n;i++) {
1274 template<
class vec_t,
class data_t>
1276 data_t &min, data_t &max) {
1283 for(
size_t i=1;i<n;i++) {
1297 template<
class vec_t,
class data_t>
1299 size_t &ix_min,
size_t &ix_max) {
1308 for(
size_t i=1;i<n;i++) {
1324 template<
class vec_t,
class data_t>
1326 size_t &ix_min, data_t &min,
1327 size_t &ix_max, data_t &max) {
1336 for(
size_t i=1;i<n;i++) {
1354 template<
class vec_t,
class data_t>
1356 size_t ix=vector_max_index<vec_t,data_t>(n,data);
1358 return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1359 }
else if (ix==n-1) {
1360 return quadratic_extremum_y<data_t>
1361 (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1363 return quadratic_extremum_y<data_t>
1364 (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1369 template<
class vec_t,
class data_t>
1371 size_t ix=vector_max_index<vec_t,data_t>(n,y);
1372 if (ix==0 || ix==n-1)
return y[ix];
1373 return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1374 y[ix-1],y[ix],y[ix+1]);
1379 template<
class vec_t,
class data_t>
1381 size_t ix=vector_max_index<vec_t,data_t>(n,y);
1382 if (ix==0 || ix==n-1)
return y[ix];
1383 return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1384 y[ix-1],y[ix],y[ix+1]);
1389 template<
class vec_t,
class data_t>
1391 size_t ix=vector_min_index<vec_t,data_t>(n,data);
1393 return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1394 }
else if (ix==n-1) {
1395 return quadratic_extremum_y<data_t>
1396 (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1398 return quadratic_extremum_y<data_t>
1399 (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1404 template<
class vec_t,
class data_t>
1406 size_t ix=vector_min_index<vec_t,data_t>(n,y);
1407 if (ix==0 || ix==n-1)
return y[ix];
1408 return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1409 y[ix-1],y[ix],y[ix+1]);
1414 template<
class vec_t,
class data_t>
1416 size_t ix=vector_min_index<vec_t,data_t>(n,y);
1417 if (ix==0 || ix==n-1)
return y[ix];
1418 return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1419 y[ix-1],y[ix],y[ix+1]);
1427 template<
class mat_t,
class data_t>
1431 std::string str=((std::string)
"Matrix with zero size (")+
1433 "matrix_max_value().";
1436 data_t max=data(0,0);
1437 for(
size_t i=0;i<m;i++) {
1438 for(
size_t j=0;j<n;j++) {
1439 if (data(i,j)>max) {
1449 template<
class mat_t,
class data_t> data_t
1451 size_t m=data.size1();
1452 size_t n=data.size2();
1454 std::string str=((std::string)
"Matrix with zero size (")+
1456 "matrix_max_value().";
1459 data_t max=data(0,0);
1460 for(
size_t i=0;i<m;i++) {
1461 for(
size_t j=0;j<n;j++) {
1462 if (data(i,j)>max) {
1472 template<
class mat_t>
double 1474 size_t m=data.size1();
1475 size_t n=data.size2();
1477 std::string str=((std::string)
"Matrix with zero size (")+
1479 "matrix_max_value_double().";
1482 double max=data(0,0);
1483 for(
size_t i=0;i<m;i++) {
1484 for(
size_t j=0;j<n;j++) {
1485 if (data(i,j)>max) {
1496 template<
class mat_t,
class data_t>
1498 size_t &i_max,
size_t &j_max, data_t &max) {
1501 std::string str=((std::string)
"Matrix with zero size (")+
1503 "matrix_max_index().";
1509 for(
size_t i=0;i<m;i++) {
1510 for(
size_t j=0;j<n;j++) {
1511 if (data(i,j)>max) {
1524 template<
class mat_t,
class data_t>
1526 size_t &i_max,
size_t &j_max, data_t &max) {
1528 size_t m=data.size1();
1529 size_t n=data.size2();
1531 std::string str=((std::string)
"Matrix with zero size (")+
1533 "matrix_max_index().";
1539 for(
size_t i=0;i<m;i++) {
1540 for(
size_t j=0;j<n;j++) {
1541 if (data(i,j)>max) {
1553 template<
class mat_t,
class data_t>
1557 std::string str=((std::string)
"Matrix with zero size (")+
1559 "matrix_min_value().";
1562 data_t min=data(0,0);
1563 for(
size_t i=0;i<m;i++) {
1564 for(
size_t j=0;j<n;j++) {
1565 if (data(i,j)<min) {
1575 template<
class mat_t,
class data_t>
1578 size_t m=data.size1();
1579 size_t n=data.size2();
1581 std::string str=((std::string)
"Matrix with zero size (")+
1583 "matrix_min_value().";
1586 data_t min=data(0,0);
1587 for(
size_t i=0;i<m;i++) {
1588 for(
size_t j=0;j<n;j++) {
1589 if (data(i,j)<min) {
1599 template<
class mat_t>
1602 size_t m=data.size1();
1603 size_t n=data.size2();
1605 std::string str=((std::string)
"Matrix with zero size (")+
1607 "matrix_min_value().";
1610 double min=data(0,0);
1611 for(
size_t i=0;i<m;i++) {
1612 for(
size_t j=0;j<n;j++) {
1613 if (data(i,j)<min) {
1624 template<
class mat_t,
class data_t>
1626 size_t &i_min,
size_t &j_min, data_t &min) {
1629 std::string str=((std::string)
"Matrix with zero size (")+
1631 "matrix_min_index().";
1637 for(
size_t i=0;i<m;i++) {
1638 for(
size_t j=0;j<n;j++) {
1639 if (data(i,j)<min) {
1652 template<
class mat_t,
class data_t>
1654 size_t &i_min,
size_t &j_min, data_t &min) {
1656 size_t m=data.size1();
1657 size_t n=data.size2();
1659 std::string str=((std::string)
"Matrix with zero size (")+
1661 "matrix_min_index().";
1667 for(
size_t i=0;i<m;i++) {
1668 for(
size_t j=0;j<n;j++) {
1669 if (data(i,j)<min) {
1681 template<
class mat_t,
class data_t>
1683 data_t &min, data_t &max) {
1690 for(
size_t i=0;i<n;i++) {
1691 for(
size_t j=0;j<m;j++) {
1692 if (data(i,j)<min) {
1694 }
else if (data(i,j)>max) {
1704 template<
class mat_t,
class data_t>
1706 data_t &min, data_t &max) {
1707 return matrix_minmax<mat_t,data_t>
1708 (data.size1(),data.size2(),data,min,max);
1714 template<
class mat_t,
class data_t>
1716 size_t &i_min,
size_t &j_min, data_t &min,
1717 size_t &i_max,
size_t &j_max, data_t &max) {
1729 for(
size_t i=0;i<n;i++) {
1730 for(
size_t j=0;j<m;j++) {
1731 if (data(i,j)<min) {
1735 }
else if (data(i,j)>max) {
1747 template<
class mat_t,
class data_t>
1751 for(
size_t i=0;i<m;i++) {
1752 for(
size_t j=0;j<n;j++) {
1761 template<
class mat_t,
class data_t>
1763 return matrix_sum<mat_t,data_t>(data.size1(),data.size2(),data);
1782 template<
class vec_t>
1785 O2SCL_ERR(
"Empty vector in function vector_lookup().",
1790 while(!std::isfinite(x[i]) && i<n-1) i++;
1796 double best=x[i], bdiff=fabs(x[i]-x0);
1798 if (std::isfinite(x[i]) && fabs(x[i]-x0)<bdiff) {
1801 bdiff=fabs(x[i]-x0);
1819 template<
class vec_t>
1829 template<
class mat_t>
1831 double x0,
size_t &i,
size_t &j) {
1833 O2SCL_ERR(
"Empty matrix in matrix_lookup().",
1837 bool found_one=
false;
1838 for(
size_t i2=0;i2<m;i2++) {
1839 for(
size_t j2=0;j2<n;j2++) {
1840 if (std::isfinite(A(i,j))) {
1841 if (found_one==
false) {
1842 dist=fabs(A(i,j)-x0);
1847 if (fabs(A(i,j)-x0)<dist) {
1848 dist=fabs(A(i,j)-x0);
1856 if (found_one==
false) {
1868 template<
class mat_t>
1870 double x0,
size_t &i,
size_t &j) {
1908 template<
class vec_t,
class data_t>
1910 size_t lo,
size_t hi) {
1912 O2SCL_ERR2(
"Low and high indexes backwards in ",
1913 "function vector_bsearch_inc().",
exc_einval);
1953 template<
class vec_t,
class data_t>
1955 size_t lo,
size_t hi) {
1957 O2SCL_ERR2(
"Low and high indexes backwards in ",
1958 "function vector_bsearch_dec().",
exc_einval);
1979 template<
class vec_t,
class data_t>
1981 size_t lo,
size_t hi) {
1982 if (x[lo]<x[hi-1]) {
1983 return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1985 return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
1995 template<
class vec_t,
class data_t>
1999 if (x[lo]<x[hi-1]) {
2000 return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
2002 return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
2021 template<
class vec_t>
2026 if (data[0]==data[1]) {
2028 }
else if (data[0]<data[1]) {
2038 for(
size_t i=0;i<n-1 && done==
false;i++) {
2039 if (data[i]!=data[i+1]) {
2053 if (data[start]<data[start+1])
return 1;
2060 if (data[start]>data[start+1]) inc=
false;
2063 for(
size_t i=start+1;i<n-1;i++) {
2065 if (data[i]>data[i+1])
return 0;
2071 for(
size_t i=start+1;i<n-1;i++) {
2072 if (data[i]<data[i+1])
return 0;
2103 template<
class vec_t>
2111 if (data[0]==data[1]) {
2113 }
else if (data[0]>data[1]) {
2118 for(
size_t i=1;i<n-1;i++) {
2120 if (data[i]>=data[i+1])
return 0;
2126 for(
size_t i=1;i<n-1;i++) {
2127 if (data[i]<=data[i+1])
return 0;
2141 template<
class vec_t>
2154 template<
class vec_t>
2156 for(
size_t i=0;i<n;i++) {
2157 if (!std::isfinite(data[i]))
return false;
2182 template<
class vec_t,
class data_t>
2186 for(
size_t i=0;i<n;i++) {
2197 template<
class vec_t,
class data_t> data_t
vector_sum(vec_t &data) {
2199 for(
size_t i=0;i<data.size();i++) {
2213 for(
size_t i=0;i<n;i++) {
2227 for(
size_t i=0;i<data.size();i++) {
2239 template<
class vec_t,
class data_t>
2247 }
else if (n == 1) {
2251 for (
size_t i = 0; i < n; i++) {
2252 const data_t xx = x[i];
2255 const data_t ax = fabs(xx);
2258 ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2261 ssq += (ax / scale) * (ax / scale);
2267 return scale * sqrt(ssq);
2273 template<
class vec_t,
class data_t> data_t
vector_norm(
const vec_t &x) {
2274 return vector_norm<vec_t,data_t>(x.size(),x);
2283 template<
class vec_t>
2291 }
else if (n == 1) {
2295 for (
size_t i = 0; i < n; i++) {
2296 const double xx = x[i];
2299 const double ax = fabs(xx);
2302 ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2305 ssq += (ax / scale) * (ax / scale);
2311 return scale * sqrt(ssq);
2317 return vector_norm_double<vec_t>(x.size(),x);
2325 template<
class vec_t,
class data_t>
2327 for(
size_t i=0;i<N;i++) {
2335 template<
class vec_t,
class data_t>
2343 template<
class mat_t,
class data_t>
2345 for(
size_t i=0;i<M;i++) {
2346 for(
size_t j=0;j<N;j++) {
2355 template<
class mat_t,
class data_t>
2366 template<
class vec_t,
class vec2_t>
2368 if (src.size()==0) {
2371 if (iout>=src.size()) {
2374 dest.resize(src.size()-1);
2376 for(
size_t i=0;i<src.size();i++) {
2398 template<
class vec_t,
class data_t>
2401 data_t *tmp=
new data_t[n];
2402 for(
size_t i=0;i<n;i++) {
2403 tmp[i]=data[(i+k)%n];
2405 for(
size_t i=0;i<n;i++) {
2417 template<
class vec_t,
class data_t>
2421 for(
size_t i=0;i<n/2;i++) {
2423 data[n-1-i]=data[i];
2434 template<
class vec_t,
class data_t>
2437 size_t n=data.size();
2439 for(
size_t i=0;i<n/2;i++) {
2441 data[n-1-i]=data[i];
2452 template<
class vec_t>
2456 for(
size_t i=0;i<n/2;i++) {
2458 data[n-1-i]=data[i];
2471 size_t n=data.size();
2473 for(
size_t i=0;i<n/2;i++) {
2475 data[n-1-i]=data[i];
2497 template<
class mat_t,
class mat_row_t>
2499 return mat_row_t(M,row);
2576 const double &operator()(
size_t row,
size_t col)
const;
2602 template<
class mat_t,
class mat_column_t>
2604 return mat_column_t(M,column);
2628 double &operator[](
size_t i) {
2629 return m_(i,column_);
2631 const double &operator[](
size_t i)
const {
2632 return m_(i,column_);
2653 m_(m), column_(column) {
2655 const double &operator[](
size_t i)
const {
2656 return m_(i,column_);
2672 template<
class vec_t>
2674 bool endline=
false) {
2679 for(
size_t i=0;i<n-1;i++) os << v[i] <<
" ";
2681 if (endline) os << std::endl;
2696 template<
class vec_t>
2697 void vector_out(std::ostream &os,
const vec_t &v,
bool endline=
false) {
2704 for(
size_t i=0;i<n-1;i++) os << v[i] <<
" ";
2706 if (endline) os << std::endl;
2712 template<
class vec_t,
class data_t>
2714 g.template vector<vec_t>(v);
2719 template<
class mat_t>
2721 for(
size_t i=0;i<M;i++) {
2722 for(
size_t j=0;j<N;j++) {
2723 if (i==j) m(i,j)=1.0;
2731 template<
class mat_t>
2746 (dat_t *v,
size_t start,
size_t last) {
2756 (
const dat_t *v,
size_t start,
size_t last) {
2772 size_t start,
size_t last) {
2775 (v,boost::numeric::ublas::range(start,last));
2790 size_t start,
size_t last) {
2793 (v,boost::numeric::ublas::range(start,last));
2809 size_t start,
size_t last) {
2812 (v,boost::numeric::ublas::range(start,last));
2825 template<
class dat_t>
2832 size_t start,
size_t last) {
2836 (v,boost::numeric::ublas::range(start,last));
2849 template<
class dat_t>
2856 size_t start,
size_t last) {
2860 (v,boost::numeric::ublas::range(start,last));
2873 template<
class dat_t>
2880 size_t start,
size_t last) {
2884 (v,boost::numeric::ublas::range(start,last));
2897 template<
class dat_t>
2904 size_t start,
size_t last) {
2908 (v,boost::numeric::ublas::range(start,last));
2935 start_(start), last_(last) {
2936 #if !O2SCL_NO_RANGE_CHECK 2938 O2SCL_ERR2(
"End before beginning in vector_range_gen::",
2939 "vector_range_gen(vec_t,size_t,size_t)",
2947 size_t last) : v_(v2.v_),
2948 start_(start+v2.start_), last_(last+v2.start_) {
2949 #if !O2SCL_NO_RANGE_CHECK 2951 O2SCL_ERR2(
"End before beginning in vector_range_gen::",
2952 "vector_range_gen(vector_range_gen,size_t,size_t)",
2955 if (last>v2.
last_) {
2956 O2SCL_ERR2(
"End beyond end of previous vector in vector_range_gen::",
2957 "vector_range_gen(vector_range_gen,size_t,size_t)",
2965 return last_-start_;
2970 #if !O2SCL_NO_RANGE_CHECK 2971 if (i+start_>=last_) {
2972 O2SCL_ERR(
"Index out of range in vector_range_gen::operator[].",
2976 return v_[i+start_];
2981 #if !O2SCL_NO_RANGE_CHECK 2982 if (i+start_>=last_) {
2987 return v_[i+start_];
3010 start_(start), last_(last) {
3011 #if !O2SCL_NO_RANGE_CHECK 3013 O2SCL_ERR2(
"End before beginning in vector_range_gen::",
3014 "vector_range_gen(vec_t,size_t,size_t)",
3022 size_t last) : v_(v2.v_),
3023 start_(start+v2.start_), last_(last+v2.start_) {
3024 #if !O2SCL_NO_RANGE_CHECK 3026 O2SCL_ERR2(
"End before beginning in vector_range_gen::",
3027 "vector_range_gen(vector_range_gen,size_t,size_t)",
3030 if (last>v2.
last_) {
3031 O2SCL_ERR2(
"End beyond end of previous vector in vector_range_gen::",
3032 "vector_range_gen(vector_range_gen,size_t,size_t)",
3040 size_t last) : v_(v2.v_),
3041 start_(start+v2.start_), last_(last+v2.start_) {
3042 #if !O2SCL_NO_RANGE_CHECK 3044 O2SCL_ERR2(
"End before beginning in vector_range_gen::",
3045 "vector_range_gen(vector_range_gen,size_t,size_t)",
3048 if (last>v2.
last_) {
3049 O2SCL_ERR2(
"End beyond end of previous vector in vector_range_gen::",
3050 "vector_range_gen(vector_range_gen,size_t,size_t)",
3058 return last_-start_;
3063 #if !O2SCL_NO_RANGE_CHECK 3064 if (i+start_>=last_) {
3069 return v_[i+start_];
3112 size_t start,
size_t last) {
3121 size_t start,
size_t last) {
3130 size_t start,
size_t last) {
3145 template<
class dat_t> std::vector<dat_t>
3147 return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3161 template<
class dat_t>
const std::vector<dat_t>
3164 return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3171 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN) 3173 #if defined (O2SCL_ARMA) || defined (DOXYGEN) 3174 #include <armadillo> 3186 template<> arma::subview_row<double>
3187 matrix_row<arma::mat,arma::subview_row<double> >
3188 (arma::mat &M,
size_t row);
3191 template<> arma::subview_col<double>
3192 matrix_column<arma::mat,arma::subview_col<double> >
3193 (arma::mat &M,
size_t column);
3200 #if defined (O2SCL_EIGEN) || defined (DOXYGEN) 3201 #include <eigen3/Eigen/Dense> 3208 double matrix_max(
const Eigen::MatrixXd &data);
3211 double matrix_min(
const Eigen::MatrixXd &data);
3214 template<> Eigen::MatrixXd::RowXpr
3215 matrix_row<Eigen::MatrixXd,Eigen::MatrixXd::RowXpr>
3216 (Eigen::MatrixXd &M,
size_t row);
3219 template<> Eigen::MatrixXd::ColXpr
3220 matrix_column<Eigen::MatrixXd,Eigen::MatrixXd::ColXpr>
3221 (Eigen::MatrixXd &M,
size_t column);
3230 #include <o2scl/vector_special.h> 3277 template<
class T,
class F,
class A>
class matrix {
data_t vector_min_quad(size_t n, const vec_t &data)
Minimum of vector by quadratic fit.
size_t vector_lookup(size_t n, const vec_t &x, double x0)
Lookup the value x0 in the first n elements of vector x.
size_t vector_min_index(size_t n, const vec_t &data)
Compute the index which holds the minimum of the first n elements of a vector.
Experimental vector range object.
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
data_t matrix_sum(size_t m, size_t n, const mat_t &data)
Compute the sum of matrix elements.
void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest)
From a given vector, create a new vector by removing a specified element.
data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector maximum by quadratic fit.
data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector minimum by quadratic fit.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void matrix_swap_rows(size_t N, mat_t &m, size_t i1, size_t i2)
Swap the first N columns of two rows in a matrix.
void matrix_make_upper(mat_t &src)
Make a matrix upper triangular by setting the lower triangular entries to zero.
bool vector_is_finite(size_t n, vec_t &data)
Test if the first n elements of a vector are finite.
Placeholder documentation of some related Boost objects.
size_t vector_bsearch_dec(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of an decreasing vector for x0.
void matrix_max_index(size_t m, size_t n, const mat_t &data, size_t &i_max, size_t &j_max, data_t &max)
Compute the maximum of a matrix and return the indices of the maximum element.
const_vector_range_gen(const vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
void vector_set_all(size_t N, vec_t &src, data_t val)
Set the first N entries in a vector to a particular value.
Generic object which represents a column of a matrix.
bool matrix_is_upper(mat_t &src)
Simple test that a matrix is upper triangular.
bool matrix_is_lower(mat_t &src)
Simple test that a matrix is lower triangular.
double & operator[](size_t i)
Return a reference ith element.
size_t vector_bsearch_inc(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of an increasing vector for x0.
A simple convenience wrapper for GSL vector objects.
void matrix_set_all(size_t M, size_t N, mat_t &src, data_t val)
Set the first (M,N) entries in a matrix to a particular value.
Generic object which represents a column of a const matrix.
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
invalid argument supplied by user
Generic object which represents a row of a matrix.
A simple matrix view object.
vector_range_gen< vec_t > vector_range(vector_range_gen< vec_t > &v, size_t start, size_t last)
Recursively create a o2scl::vector_range_gen object from a vector range.
void vector_reverse(size_t n, vec_t &data)
Reverse the first n elements of a vector.
size_t last_
The end() iterator.
void vector_rotate(size_t n, vec_t &data, size_t k)
"Rotate" a vector so that the kth element is now the beginning
void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2)
Swap of of the first N elements of two double-precision vectors.
size_t last_
The end() iterator.
double matrix_min_value_double(const mat_t &data)
Compute the minimum of a matrix.
The default matrix type from uBlas.
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
Generic object which represents a row of a matrix.
const_matrix_row_gen(const mat_t &m, size_t row)
Create a row object from row row of matrix m.
double vector_norm_double(size_t n, const vec_t &x)
Compute the norm of the first n entries of a vector of double precision numbers.
void matrix_make_lower(mat_t &src)
Make a matrix lower triangular by setting the upper triangular entries to zero.
double matrix_max_value_double(const mat_t &data)
Compute the maximum of a matrix.
data_t vector_sum(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector.
void vector_swap(size_t N, vec_t &v1, vec2_t &v2)
Swap the first N elements of two vectors.
void vector_largest(size_t n, vec_t &data, size_t k, vec_t &largest)
Find the k largest entries of the first n elements of a vector.
std::vector< dat_t > vector_range_copy(const std::vector< dat_t > &v, size_t start, size_t last)
Vector range function template for std::vector
const_vector_range_gen(const vector_range_gen< vec_t > &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
matrix_row_gen(mat_t &m, size_t row)
Create a row object from row row of matrix m.
vec_t & v_
A reference to the original vector.
void matrix_swap_double(size_t M, size_t N, mat_t &m1, mat2_t &m2)
Swap of the first elements in two double-precision matrices.
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts the first n elements of a vector (in increasing order) ...
void matrix_swap_cols(size_t M, mat_t &m, size_t j1, size_t j2)
Swap the first M rows of two columns in a matrix.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
size_t size() const
Return the vector size.
void matrix_minmax_index(size_t n, size_t m, const mat_t &data, size_t &i_min, size_t &j_min, data_t &min, size_t &i_max, size_t &j_max, data_t &max)
Compute the minimum and maximum of a matrix and return their locations.
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
const vec_t & v_
A reference to the original vector.
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
const_vector_range_gen(const const_vector_range_gen &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
void vector_minmax_value(size_t n, vec_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
void matrix_set_identity(size_t M, size_t N, mat_t &m)
Set a matrix to unity on the diagonal and zero otherwise.
data_t vector_max_quad(size_t n, const vec_t &data)
Maximum of vector by quadratic fit.
vector_range_gen(const vector_range_gen &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
mat_column_t matrix_column(mat_t &M, size_t column)
Construct a column of a matrix.
size_t start_
The index offset.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
void vector_smallest(size_t n, vec_t &data, size_t k, vec_t &smallest)
Find the k smallest entries of the first n elements of a vector.
void vector_grid(uniform_grid< data_t > g, vec_t &v)
Fill a vector with a specified grid.
void vector_minmax(size_t n, vec_t &data, size_t &ix_min, data_t &min, size_t &ix_max, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
data_t vector_norm(size_t n, const vec_t &x)
Compute the norm of the first n entries of a vector of floating-point (single or double precision) nu...
A simple convenience wrapper for GSL matrix objects.
size_t vector_max_index(size_t n, const vec_t &data)
Compute the index which holds the maximum of the first n elements of a vector.
void vector_max(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the maximum of the first n elements of a vector.
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
void vector_minmax_index(size_t n, vec_t &data, size_t &ix_min, size_t &ix_max)
Compute the minimum and maximum of the first n elements of a vector.
void vector_smallest_index(size_t n, const vec_t &data, size_t k, vec_size_t &index)
Find the indexes of the k smallest entries among the first n entries of a vector. ...
double matrix_max(const arma::mat &data)
Armadillo version of matrix_max()
const double & operator[](size_t i) const
Return a const reference ith element.
const double & operator[](size_t i) const
Return a const reference ith element.
double vector_sum_double(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector of double-precision numbers.
vector_range_gen(vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
void matrix_minmax(size_t n, size_t m, const mat_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of a matrix.
Experimental const vector range object.
size_t vector_bsearch(const data_t x0, const vec_t &x, size_t lo, size_t hi)
Binary search a part of a monotonic vector for x0.
void matrix_transpose(mat_t &src, mat2_t &dest)
Simple transpose.
void sort_index_downheap(size_t N, const vec_t &data, vec_size_t &order, size_t k)
Provide a downheap() function for vector_sort_index()
void vector_reverse_double(size_t n, vec_t &data)
Reverse the first n elements in a vector of double precision numbers.
int vector_is_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are monotonic and increasing or decreasing.
int vector_is_strictly_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are strictly monotonic and determine if they are increasing ...
void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2)
Swap the first M rows of two columns in a double-precision matrix.
void matrix_min_index(size_t n, size_t m, const mat_t &data, size_t &i_min, size_t &j_min, data_t &min)
Compute the minimum of a matrix and return the indices of the minimum element.
double matrix_min(const arma::mat &data)
Armadillo version of matrix_min()
void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2)
Swap of the first elements in two matrices.
void matrix_swap_rows_double(size_t N, mat_t &m, size_t i1, size_t i2)
Swap the first N columns of two rows in a double-precision matrix.
data_t matrix_max_value(size_t m, const size_t n, const mat_t &data)
Compute the maximum of the lower-left part of a matrix.
size_t start_
The index offset.
std::string itos(int x)
Convert an integer to a string.
std::string szttos(size_t x)
Convert a size_t to a string.
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
void matrix_lookup(size_t m, size_t n, const mat_t &A, double x0, size_t &i, size_t &j)
Lookup an element in the first $(m,n)$ entries in a matrix.
double & operator[](size_t i)
Return a reference to the ith column of the selected row.
data_t matrix_min_value(size_t m, size_t n, const mat_t &data)
Compute the minimum of a matrix.
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
The default vector type from uBlas.
size_t size() const
Return the vector size.
void vector_sort(size_t n, vec_t &data)
Sort a vector (in increasing order)
void sort_downheap(vec_t &data, size_t n, size_t k)
Provide a downheap() function for vector_sort()