vector.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_VECTOR_H
24 #define O2SCL_VECTOR_H
25 
26 /** \file vector.h
27  \brief Assorted generic vector functions
28 
29  This file contains a set of template functions which can be
30  applied to almost any vector or matrix type which allow element
31  access through <tt>operator[]</tt> (for vectors) or
32  <tt>operator(,)</tt> for matrices. Detailed requirements on the
33  template parameters are given in the functions below.
34 
35  For a general discussion of vectors and matrices in \o2, see the
36  \ref vecmat_section of the User's Guide.
37 
38  For statistics operations not included here, see \ref vec_stats.h
39  in the directory \c src/other . Also related are the matrix output
40  functions, \ref o2scl::matrix_out(), which is defined in \ref
41  columnify.h because they utilize the class \ref o2scl::columnify to
42  format the output.
43 
44  For functions which search for a value in an ordered (either
45  increasing or decreasing) vector, see the class \ref
46  o2scl::search_vec .
47 */
48 #include <iostream>
49 #include <cmath>
50 #include <string>
51 #include <fstream>
52 #include <sstream>
53 
54 #include <gsl/gsl_vector.h>
55 #include <gsl/gsl_sys.h>
56 #include <gsl/gsl_matrix.h>
57 
58 #include <o2scl/misc.h>
59 #include <o2scl/uniform_grid.h>
60 
61 namespace o2scl {
62 
63  /** \brief A simple convenience wrapper for GSL vector objects
64 
65  \warning This uses typecasts on externally allocated GSL
66  pointers and is not safe or fully const-correct.
67  */
69  const double *d;
70  size_t sz;
71  public:
72  gsl_vector_wrap(gsl_vector *m) {
73  d=(const double *)m->data;
74  sz=m->size;
75  }
76  double operator[](size_t i) const {
77  return d[i];
78  }
79  size_t size() {
80  return sz;
81  }
82  };
83 
84  /** \brief A simple convenience wrapper for GSL matrix objects
85 
86  \warning This uses typecasts on externally allocated GSL
87  pointers and is not safe or fully const-correct.
88  */
90  const double *d;
91  size_t sz1;
92  size_t sz2;
93  size_t tda;
94  public:
95  gsl_matrix_wrap(gsl_matrix *m) {
96  d=(const double *)m->data;
97  sz1=m->size1;
98  sz2=m->size2;
99  tda=m->tda;
100  }
101  double operator()(size_t i, size_t j) const {
102  return *(d+i*tda+j);
103  }
104  size_t size1() const {
105  return sz1;
106  }
107  size_t size2() const {
108  return sz2;
109  }
110  };
111 
112  /// \name Copying vectors and matrices in src/base/vector.h
113  //@{
114  /** \brief Simple vector copy
115 
116  Copy \c src to \c dest, resizing \c dest if it is too small
117  to hold <tt>src.size()</tt> elements.
118 
119  This function will work for any classes \c vec_t and
120  \c vec2_t which have suitably defined <tt>operator[]</tt>,
121  <tt>size()</tt>, and <tt>resize()</tt> methods.
122  */
123  template<class vec_t, class vec2_t>
124  void vector_copy(const vec_t &src, vec2_t &dest) {
125  size_t N=src.size();
126  if (dest.size()<N) dest.resize(N);
127  size_t i, m=N%4;
128  for(i=0;i<m;i++) {
129  dest[i]=src[i];
130  }
131  for(i=m;i+3<N;i+=4) {
132  dest[i]=src[i];
133  dest[i+1]=src[i+1];
134  dest[i+2]=src[i+2];
135  dest[i+3]=src[i+3];
136  }
137  return;
138  }
139 
140  /** \brief Simple vector copy of the first N elements
141 
142  Copy the first \c N elements of \c src to \c dest.
143  It is assumed that the memory allocation for \c dest
144  has already been performed.
145 
146  This function will work for any class <tt>vec2_t</tt> which has
147  an operator[] which returns a reference to the corresponding
148  element and class <tt>vec_t</tt> with an operator[] which
149  returns either a reference or the value of the corresponding
150  element.
151  */
152  template<class vec_t, class vec2_t>
153  void vector_copy(size_t N, const vec_t &src, vec2_t &dest) {
154  size_t i, m=N%4;
155  for(i=0;i<m;i++) {
156  dest[i]=src[i];
157  }
158  for(i=m;i+3<N;i+=4) {
159  dest[i]=src[i];
160  dest[i+1]=src[i+1];
161  dest[i+2]=src[i+2];
162  dest[i+3]=src[i+3];
163  }
164  return;
165  }
166 
167  /** \brief Simple matrix copy
168 
169  Copy \c src to \c dest, resizing \c dest if it is too small.
170 
171  This function will work for any classes \c mat_t and
172  \c mat2_t which have suitably defined <tt>operator()</tt>,
173  <tt>size()</tt>, and <tt>resize()</tt> methods.
174  */
175  template<class mat_t, class mat2_t>
176  void matrix_copy(mat_t &src, mat2_t &dest) {
177  size_t m=src.size1();
178  size_t n=src.size2();
179  if (dest.size1()<m || dest.size2()<n) dest.resize(m,n);
180  for(size_t i=0;i<m;i++) {
181  for(size_t j=0;j<n;j++) {
182  dest(i,j)=src(i,j);
183  }
184  }
185  }
186 
187  /** \brief Simple matrix copy of the first \f$ (M,N) \f$
188  matrix elements
189 
190  Copy the first <tt>(M,N)</tt> elements of \c src to \c dest. It
191  is assumed that the memory allocation for \c dest has already
192  been performed.
193 
194  This function will work for any class <tt>vec2_t</tt> which has
195  an operator[][] which returns a reference to the corresponding
196  element and class <tt>vec_t</tt> with an operator[][] which
197  returns either a reference or the value of the corresponding
198  element.
199  */
200  template<class mat_t, class mat2_t>
201  void matrix_copy(size_t M, size_t N, mat_t &src, mat2_t &dest) {
202  for(size_t i=0;i<M;i++) {
203  for(size_t j=0;j<N;j++) {
204  dest(i,j)=src(i,j);
205  }
206  }
207  }
208  //@}
209 
210  /// \name Tranpositions in src/base/vector.h
211  //@{
212  /** \brief Simple transpose
213 
214  Copy the transpose of \c src to \c dest, resizing \c dest if it
215  is too small.
216 
217  This function will work for any classes \c mat_t and
218  \c mat2_t which have suitably defined <tt>operator()</tt>,
219  <tt>size()</tt>, and <tt>resize()</tt> methods.
220  */
221  template<class mat_t, class mat2_t>
222  void matrix_transpose(mat_t &src, mat2_t &dest) {
223  size_t m=src.size1();
224  size_t n=src.size2();
225  if (dest.size1()<n || dest.size2()<m) dest.resize(n,m);
226  for(size_t i=0;i<m;i++) {
227  for(size_t j=0;j<n;j++) {
228  dest(i,j)=src(j,i);
229  }
230  }
231  }
232 
233  /** \brief Simple transpose of the first \f$ (m,n) \f$
234  matrix elements
235 
236  Copy the transpose of the first \c m rows and the first \c cols
237  of the matrix \c src into the matrix \c dest
238 
239  This function will work for any classes \c mat_t and \c mat2_t
240  which has a suitably defined <tt>operator()</tt> method.
241  */
242  template<class mat_t, class mat2_t>
243  void matrix_transpose(size_t m, size_t n, mat_t &src, mat2_t &dest) {
244  for(size_t i=0;i<m;i++) {
245  for(size_t j=0;j<n;j++) {
246  dest(i,j)=src(j,i);
247  }
248  }
249  }
250 
251  /** \brief Simple in-place transpose
252 
253  Transpose the matrix \c src . If the matrix is not square,
254  only the upper-left square part of the matrix will be
255  transposed.
256 
257  This function will work for any classes \c mat_t and
258  \c mat2_t which have suitably defined <tt>operator()</tt>,
259  <tt>size()</tt>, and <tt>resize()</tt> methods.
260  */
261  template<class mat_t, class data_t>
262  void matrix_transpose(mat_t &src) {
263  size_t m=src.size1();
264  size_t n=src.size2();
265  // Choose the smaller of n and m
266  if (m<n) n=m;
267  data_t tmp;
268  for(size_t i=0;i<n;i++) {
269  for(size_t j=0;j<n;j++) {
270  tmp=src(i,j);
271  src(i,j)=src(j,i);
272  src(j,i)=tmp;
273  }
274  }
275  }
276 
277  /** \brief Simple in-place transpose of the first \f$ (m,n) \f$
278  matrix elements
279 
280  Copy the transpose of the first \c m rows and the first \c cols
281  of the matrix \c src into the matrix \c dest
282 
283  This function will work for any classes \c mat_t and \c mat2_t
284  which has a suitably defined <tt>operator()</tt> method.
285  */
286  template<class mat_t, class data_t>
287  void matrix_transpose(size_t m, size_t n, mat_t &src) {
288  data_t tmp;
289  for(size_t i=0;i<m;i++) {
290  for(size_t j=0;j<n;j++) {
291  tmp=src(i,j);
292  src(i,j)=src(j,i);
293  src(j,i)=tmp;
294  }
295  }
296  }
297  //@}
298 
299  /// \name Upper and lower triangular functions in src/base/vector.h
300  //@{
301  /** \brief Simple test that a matrix is lower triangular
302  */
303  template<class mat_t> bool matrix_is_lower(mat_t &src) {
304  size_t m=src.size1();
305  size_t n=src.size2();
306  bool ret=true;
307  for(size_t i=0;i<m;i++) {
308  for(size_t j=i+1;j<n;j++) {
309  if (src(i,j)!=0.0) return false;
310  }
311  }
312  return ret;
313  }
314 
315  /** \brief Simple test that a matrix is upper triangular
316  */
317  template<class mat_t> bool matrix_is_upper(mat_t &src) {
318  size_t m=src.size1();
319  size_t n=src.size2();
320  bool ret=true;
321  for(size_t j=0;j<n;j++) {
322  for(size_t i=j+1;i<m;i++) {
323  if (src(i,j)!=0.0) return false;
324  }
325  }
326  return ret;
327  }
328 
329  /** \brief Make a matrix lower triangular by setting the upper
330  triangular entries to zero
331  */
332  template<class mat_t> void matrix_make_lower(mat_t &src) {
333  size_t m=src.size1();
334  size_t n=src.size2();
335  for(size_t i=0;i<m;i++) {
336  for(size_t j=i+1;j<n;j++) {
337  src(i,j)=0.0;
338  }
339  }
340  return;
341  }
342 
343  /** \brief Make a matrix upper triangular by setting the lower
344  triangular entries to zero
345  */
346  template<class mat_t> void matrix_make_upper(mat_t &src) {
347  size_t m=src.size1();
348  size_t n=src.size2();
349  for(size_t j=0;j<n;j++) {
350  for(size_t i=j+1;i<m;i++) {
351  src(i,j)=0.0;
352  }
353  }
354  return;
355  }
356 
357  /** \brief Simple test that a matrix is lower triangular
358  for the first \c m rows and \c n columns
359  */
360  template<class mat_t> bool matrix_is_lower(size_t m, size_t n,
361  mat_t &src) {
362  bool ret=true;
363  for(size_t i=0;i<m;i++) {
364  for(size_t j=i+1;j<n;j++) {
365  if (src(i,j)!=0.0) return false;
366  }
367  }
368  return ret;
369  }
370 
371  /** \brief Simple test that a matrix is upper triangular
372  for the first \c m rows and \c n columns
373  */
374  template<class mat_t> bool matrix_is_upper(size_t m, size_t n,
375  mat_t &src) {
376  bool ret=true;
377  for(size_t j=0;j<n;j++) {
378  for(size_t i=j+1;i<m;i++) {
379  if (src(i,j)!=0.0) return false;
380  }
381  }
382  return ret;
383  }
384 
385  /** \brief Make the first \c m rows and \c n columns of a matrix
386  lower triangular by setting the upper triangular entries to zero
387  */
388  template<class mat_t> void matrix_make_lower(size_t m, size_t n,
389  mat_t &src) {
390  for(size_t i=0;i<m;i++) {
391  for(size_t j=i+1;j<n;j++) {
392  src(i,j)=0.0;
393  }
394  }
395  return;
396  }
397 
398  /** \brief Make the first \c m rows and \c n columns of a matrix
399  upper triangular by setting the lower triangular entries to zero
400  */
401  template<class mat_t> void matrix_make_upper(size_t m, size_t n,
402  mat_t &src) {
403  for(size_t j=0;j<n;j++) {
404  for(size_t i=j+1;i<m;i++) {
405  src(i,j)=0.0;
406  }
407  }
408  return;
409  }
410  //@}
411 
412  /// \name Swapping parts of vectors and matrices in src/base/vector.h
413  //@{
414  /** \brief Swap the first \c N elements of two vectors
415 
416  This function swaps the elements of \c v1 and \c v2, one element
417  at a time.
418  */
419  template<class vec_t, class vec2_t, class data_t>
420  void vector_swap(size_t N, vec_t &v1, vec2_t &v2) {
421  data_t temp;
422  size_t i, m=N%4;
423  for(i=0;i<m;i++) {
424  temp=v1[i];
425  v1[i]=v2[i];
426  v2[i]=temp;
427  }
428  for(i=m;i+3<N;i+=4) {
429  temp=v1[i];
430  v1[i]=v2[i];
431  v2[i]=temp;
432  temp=v1[i+1];
433  v1[i+1]=v2[i+1];
434  v2[i+1]=temp;
435  temp=v1[i+2];
436  v1[i+2]=v2[i+2];
437  v2[i+2]=temp;
438  temp=v1[i+3];
439  v1[i+3]=v2[i+3];
440  v2[i+3]=temp;
441  }
442  return;
443  }
444 
445  /** \brief Swap all elements in two vectors
446 
447  This function swaps the elements of \c v1 and \c v2, one element
448  at a time.
449 
450  \note It is almost always better to use <tt>std::swap</tt>
451  than this function, which is provided only in cases where
452  one knows one is going to be forced to use a vector type
453  without a properly defined <tt>std::swap</tt> method.
454  */
455  template<class vec_t, class vec2_t, class data_t>
456  void vector_swap(vec_t &v1, vec2_t &v2) {
457  size_t N=v1.size();
458  if (v2.size()<N) N=v2.size();
459  data_t temp;
460  size_t i, m=N%4;
461  for(i=0;i<m;i++) {
462  temp=v1[i];
463  v1[i]=v2[i];
464  v2[i]=temp;
465  }
466  for(i=m;i+3<N;i+=4) {
467  temp=v1[i];
468  v1[i]=v2[i];
469  v2[i]=temp;
470  temp=v1[i+1];
471  v1[i+1]=v2[i+1];
472  v2[i+1]=temp;
473  temp=v1[i+2];
474  v1[i+2]=v2[i+2];
475  v2[i+2]=temp;
476  temp=v1[i+3];
477  v1[i+3]=v2[i+3];
478  v2[i+3]=temp;
479  }
480  return;
481  }
482 
483  /** \brief Swap of of the first N elements of two
484  double-precision vectors
485 
486  This function swaps the elements of \c v1 and \c v2, one element
487  at a time.
488  */
489  template<class vec_t, class vec2_t>
490  void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2) {
491  return vector_swap<vec_t,vec2_t,double>(N,v1,v2);
492  }
493 
494  /** \brief Swap of all the elements in two
495  double-precision vectors
496 
497  This function swaps the elements of \c v1 and \c v2, one element
498  at a time.
499 
500  \note It is almost always better to use <tt>std::swap</tt>
501  than this function, which is provided only in cases where
502  one knows one is going to be forced to use a vector type
503  without a properly defined <tt>std::swap</tt> method.
504  */
505  template<class vec_t, class vec2_t>
506  void vector_swap_double(vec_t &v1, vec2_t &v2) {
507  return vector_swap<vec_t,vec2_t,double>(v1,v2);
508  }
509 
510  /** \brief Swap two elements in a vector
511 
512  This function swaps the element \c i and element \c j of vector
513  \c v1.
514  */
515  template<class vec_t, class data_t>
516  void vector_swap(vec_t &v, size_t i, size_t j) {
517  data_t temp=v[i];
518  v[i]=v[j];
519  v[j]=temp;
520  return;
521  }
522 
523  /** \brief Swap two elements in a double-precision vector
524 
525  This function swaps the element \c i and element \c j of vector
526  \c v1.
527 
528  This function is used in \ref o2scl_linalg::QRPT_decomp() .
529  */
530  template<class vec_t>
531  void vector_swap_double(vec_t &v, size_t i, size_t j) {
532  return vector_swap<vec_t,double>(v,i,j);
533  }
534 
535  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
536  matrices
537 
538  This function swaps the elements of \c v1 and \c v2, one element
539  at a time.
540  */
541  template<class mat_t, class mat2_t, class data_t>
542  void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2) {
543  data_t temp;
544  for(size_t i=0;i<M;i++) {
545  for(size_t j=0;j<N;j++) {
546  temp=v1[i][j];
547  v1[i][j]=v2[i][j];
548  v2[i][j]=temp;
549  }
550  }
551  return;
552  }
553 
554  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
555  double-precision matrices
556 
557  This function swaps the elements of \c m1 and \c m2, one element
558  at a time.
559  */
560  template<class mat_t, class mat2_t, class data_t>
561  void matrix_swap_double(size_t M, size_t N, mat_t &m1, mat2_t &m2) {
562  return matrix_swap<mat_t,mat2_t,double>(M,N,m1,m2);
563  }
564 
565  /** \brief Swap two elements in a matrix
566 
567  This function swaps the element <tt>(i1,j1)</tt> and
568  element <tt>(i2,j2)</tt> of matrix \c m1.
569  */
570  template<class mat_t, class data_t>
571  void matrix_swap(mat_t &m, size_t i1, size_t j1, size_t i2, size_t j2) {
572  data_t temp=m(i1,j1);
573  m(i1,j1)=m(i2,j2);
574  m(i2,j2)=temp;
575  return;
576  }
577 
578  /** \brief Swap two elements in a double-precision matrix
579 
580  This function swaps the element \c i and element \c j of matrix
581  \c v1.
582  */
583  template<class mat_t>
584  void matrix_swap_double(mat_t &m, size_t i1, size_t j1,
585  size_t i2, size_t j2) {
586  return matrix_swap<mat_t,double>(m,i1,j1,i2,j2);
587  }
588 
589  /** \brief Swap the first \c M rows of two columns in a matrix
590 
591  This function swaps the element <tt>(i1,j1)</tt> and
592  element <tt>(i2,j2)</tt> of matrix \c m1.
593  */
594  template<class mat_t, class data_t>
595  void matrix_swap_cols(size_t M, mat_t &m, size_t j1, size_t j2) {
596  data_t temp;
597  for(size_t i=0;i<M;i++) {
598  temp=m(i,j1);
599  m(i,j1)=m(i,j2);
600  m(i,j2)=temp;
601  }
602  return;
603  }
604 
605  /** \brief Swap the first \c M rows of two columns in a
606  double-precision matrix
607 
608  This function swaps the element \c i and element \c j of matrix
609  \c v1.
610  */
611  template<class mat_t>
612  void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2) {
613  return matrix_swap_cols<mat_t,double>(M,m,j1,j2);
614  }
615 
616  /** \brief Swap the first \c N columns of two rows in a
617  matrix
618 
619  This function swaps the element <tt>(i1,j1)</tt> and
620  element <tt>(i2,j2)</tt> of matrix \c m1.
621  */
622  template<class mat_t, class data_t>
623  void matrix_swap_rows(size_t N, mat_t &m, size_t i1, size_t i2) {
624  data_t temp;
625  for(size_t j=0;j<N;j++) {
626  temp=m(i1,j);
627  m(i1,j)=m(i2,j);
628  m(i2,j)=temp;
629  }
630  return;
631  }
632 
633  /** \brief Swap the first \c N columns of two rows in a
634  double-precision matrix
635 
636  This function swaps the element \c i and element \c j of matrix
637  \c v1.
638  */
639  template<class mat_t>
640  void matrix_swap_rows_double(size_t N, mat_t &m, size_t i1, size_t i2) {
641  return matrix_swap_rows<mat_t,double>(N,m,i1,i2);
642  }
643  //@}
644 
645  /// \name Sorting vectors in src/base/vector.h
646  //@{
647  /** \brief Provide a downheap() function for vector_sort()
648  */
649  template<class vec_t, class data_t>
650  void sort_downheap(vec_t &data, size_t n, size_t k) {
651 
652  data_t v=data[k];
653 
654  while (k<=n/2) {
655  size_t j=2*k;
656 
657  if (j<n && data[j] < data[j+1]) j++;
658  if (!(v < data[j])) break;
659  data[k]=data[j];
660  k=j;
661  }
662 
663  data[k]=v;
664 
665  return;
666  }
667 
668  /** \brief Sort a vector (in increasing order)
669 
670  This is a generic sorting template function using a heapsort
671  algorithm. It will work for any types \c data_t and \c vec_t for
672  which
673  - \c data_t has a non-const version of <tt>operator=</tt>
674  - \c data_t has a less than operator to compare elements
675  - <tt>vec_t::operator[]</tt> returns a non-const reference
676  to an object of type \c data_t
677 
678  In particular, it will work with the STL template class
679  <tt>std::vector</tt>, and arrays and pointers of numeric,
680  character, and string objects.
681 
682  For example,
683  \code
684  std::string list[3]={"dog","cat","fox"};
685  vector_sort<std::string[3],std::string>(3,list);
686  \endcode
687 
688  \note With this function template alone, the user cannot avoid
689  explicitly specifying the template types for this function
690  because there is no parameter of type \c data_t, and function
691  templates cannot handle default template types. For this
692  reason, the function template \ref o2scl::vector_sort_double() was
693  also created which provides the convenience of not requiring
694  the user to specify the vector template type.
695 
696  \note This sorting routine is not stable, i.e. equal elements
697  have arbtrary final ordering
698 
699  \note If \c n is zero, this function will do nothing and will
700  not call the error handler.
701 
702  This works similarly to the GSL function <tt>gsl_sort_vector()</tt>.
703  */
704  template<class vec_t, class data_t>
705  void vector_sort(size_t n, vec_t &data) {
706 
707  size_t N;
708  size_t k;
709 
710  if (n==0) return;
711 
712  N=n-1;
713  k=N/2;
714  k++;
715  do {
716  k--;
717  sort_downheap<vec_t,data_t>(data,N,k);
718  } while (k > 0);
719 
720  while (N > 0) {
721  data_t tmp=data[0];
722  data[0]=data[N];
723  data[N]=tmp;
724  N--;
725  sort_downheap<vec_t,data_t>(data,N,0);
726  }
727 
728  return;
729  }
730 
731  /** \brief Provide a downheap() function for vector_sort_index()
732  */
733  template<class vec_t, class vec_size_t>
734  void sort_index_downheap(size_t N, const vec_t &data, vec_size_t &order,
735  size_t k) {
736 
737  const size_t pki = order[k];
738 
739  while (k <= N / 2) {
740  size_t j = 2 * k;
741 
742  if (j < N && data[order[j]] < data[order[j + 1]]) {
743  j++;
744  }
745 
746  // [GSL] Avoid infinite loop if nan
747  if (!(data[pki] < data[order[j]])) {
748  break;
749  }
750 
751  order[k] = order[j];
752 
753  k = j;
754  }
755 
756  order[k] = pki;
757 
758  return;
759  }
760 
761  /** \brief Create a permutation which sorts
762  the first \c n elements of a vector (in increasing order)
763 
764  This function takes a vector \c data and arranges a list of
765  indices in \c order, which give a sorted version of the vector.
766  The value <tt>order[i]</tt> gives the index of entry in in \c
767  data which corresponds to the <tt>i</tt>th value in the sorted
768  vector. The vector \c data is unchanged by this function, and
769  the initial values in \c order are ignored. Before calling this
770  function, \c order must already be allocated as a vector of size
771  \c n.
772 
773  For example, after calling this function, a sorted version the
774  vector can be output with
775  \code
776  size_t n=5;
777  double data[5]={3.1,4.1,5.9,2.6,3.5};
778  permutation order(n);
779  vector_sort_index(n,data,order);
780  for(size_t i=0;i<n;i++) {
781  cout << data[order[i]] << endl;
782  }
783  \endcode
784 
785  To create a permutation which stores as its <tt>i</tt>th element,
786  the index of <tt>data[i]</tt> in the sorted vector, you can
787  invert the permutation created by this function.
788 
789  This is a generic sorting template function. It will work for
790  any types \c vec_t and \c vec_size_t for which
791  - \c vec_t has an <tt>operator[]</tt>, and
792  - \c vec_size_t has an <tt>operator[]</tt> which returns
793  a \c size_t .
794  One possible type for \c vec_size_t is \ref o2scl::permutation.
795 
796  This works similarly to the GSL function <tt>gsl_sort_index()</tt>.
797 
798  */
799  template<class vec_t, class vec_size_t>
800  void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order) {
801  size_t N;
802  size_t i, k;
803 
804  if (n == 0) return;
805 
806  // [GSL] Set permutation to identity
807 
808  for (i = 0 ; i < n ; i++) {
809  order[i] = i;
810  }
811 
812  /* [GSL] We have n_data elements, last element is at 'n_data-1',
813  first at '0' Set N to the last element number.
814  */
815  N = n - 1;
816 
817  k = N / 2;
818  // [GSL] Compensate the first use of 'k--'
819  k++;
820  do {
821  k--;
822  sort_index_downheap<vec_t,vec_size_t>(N,data,order,k);
823  } while (k > 0);
824 
825  while (N > 0) {
826 
827  // [GSL] First swap the elements
828  size_t tmp = order[0];
829  order[0] = order[N];
830  order[N] = tmp;
831 
832  // [GSL] Then process the heap
833  N--;
834 
835  sort_index_downheap<vec_t,vec_size_t>(N,data,order,0);
836  }
837 
838  return;
839  }
840 
841  /** \brief Create a permutation which sorts a vector
842  (in increasing order)
843 
844  This function takes a vector \c data and arranges a list of
845  indices in \c order, which give a sorted version of the vector.
846  The value <tt>order[i]</tt> gives the index of entry in in \c
847  data which corresponds to the <tt>i</tt>th value in the sorted
848  vector. The vector \c data is unchanged by this function, and
849  the initial values in \c order are ignored. Before calling this
850  function, \c order must already be allocated as a vector of size
851  \c n.
852 
853  For example, after calling this function, a sorted version the
854  vector can be output with
855  \code
856  size_t n=5;
857  double data[5]={3.1,4.1,5.9,2.6,3.5};
858  permutation order(n);
859  vector_sort_index(n,data,order);
860  for(size_t i=0;i<n;i++) {
861  cout << data[order[i]] << endl;
862  }
863  \endcode
864 
865  To create a permutation which stores as its <tt>i</tt>th element,
866  the index of <tt>data[i]</tt> in the sorted vector, you can
867  invert the permutation created by this function.
868 
869  This is a generic sorting template function. It will work for
870  any types \c vec_t and \c vec_size_t for which
871  - \c vec_t has an <tt>operator[]</tt>, and
872  - \c vec_size_t has an <tt>operator[]</tt> which returns
873  a \c size_t .
874  One possible type for \c vec_size_t is \ref o2scl::permutation.
875 
876  This works similarly to the GSL function <tt>gsl_sort_index()</tt>.
877 
878  */
879  template<class vec_t, class vec_size_t>
880  void vector_sort_index(const vec_t &data, vec_size_t &order) {
881  vector_sort_index(data.size(),data,order);
882  return;
883  }
884 
885  /** \brief Sort a vector of doubles (in increasing order)
886 
887  This function is just a wrapper for
888  \code
889  vector_sort<vec_t,double>(n,data);
890  \endcode
891  See the documentation of \ref o2scl::vector_sort() for more
892  details.
893  */
894  template<class vec_t>
895  void vector_sort_double(size_t n, vec_t &data) {
896  return vector_sort<vec_t,double>(n,data);
897  }
898  //@}
899 
900  /// \name Smallest or largest subset functions in src/base/vector.h
901  //@{
902  /** \brief Find the k smallest entries of the first \c n elements
903  of a vector
904 
905  Given a vector \c data of size \c n, this function sets the
906  first \c k entries of the vector \c smallest to the k smallest
907  entries from vector \c data in ascending order. The vector \c
908  smallest must be allocated beforehand to hold at least \c k
909  elements.
910 
911  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
912 
913  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
914  \f$ k << N \f$.
915 
916  If \c k is zero, then this function does nothing and
917  returns \ref o2scl::success .
918  */
919  template<class vec_t, class data_t>
920  void vector_smallest(size_t n, vec_t &data, size_t k, vec_t &smallest) {
921  if (k>n) {
922  O2SCL_ERR2("Subset length greater than size in ",
923  "function vector_smallest().",exc_einval);
924  }
925  if (k==0 || n==0) {
926  O2SCL_ERR2("Vector size zero or k zero in ",
927  "function vector_smallest().",exc_einval);
928  }
929 
930  // Take the first element
931  size_t j=1;
932  data_t xbound=data[0];
933  smallest[0]=xbound;
934 
935  // Examine the remaining elements
936  for(size_t i=1;i<n;i++) {
937  data_t xi=data[i];
938  if (j<k) {
939  j++;
940  } else if (xi>=xbound) {
941  continue;
942  }
943  size_t i1;
944  for(i1=j-1;i1>0;i1--) {
945  if (xi>smallest[i1-1]) break;
946  smallest[i1]=smallest[i1-1];
947  }
948  smallest[i1]=xi;
949  xbound=smallest[j-1];
950  }
951  return;
952  }
953 
954  /** \brief Find the k smallest entries of a vector
955  of a vector
956 
957  Given a vector \c data, this function sets the first \c k
958  entries of the vector \c smallest to the k smallest entries from
959  vector \c data in ascending order. The vector \c smallest
960  is resized if necessary to hold at least \c k elements.
961 
962  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
963 
964  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
965  \f$ k << N \f$.
966 
967  If \c k is zero, then this function does nothing and
968  returns \ref o2scl::success .
969  */
970  template<class vec_t, class data_t>
971  void vector_smallest(vec_t &data, size_t k, vec_t &smallest) {
972  size_t n=data.size();
973  if (smallest.size()<k) smallest.resize(k);
974  return vector_smallest(n,data,k,smallest);
975  }
976 
977  /** \brief Find the indexes of the k smallest entries among the
978  first \c n entries of a vector
979 
980  Given a vector \c data, this function sets the first \c k
981  entries of the vector \c smallest equal to the indexes of the k
982  smallest entries from vector \c data in ascending order. The
983  vector \c smallest is resized if necessary to hold at least \c k
984  elements.
985 
986  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
987  \f$ k << N \f$.
988 
989  If \c k is zero or \c n is zero or \f$ k > n\f$, then this
990  function calls the error handler.
991  */
992  template<class vec_t, class data_t, class vec_size_t>
993  void vector_smallest_index(size_t n, const vec_t &data, size_t k,
994  vec_size_t &index) {
995  if (k>n) {
996  O2SCL_ERR2("Subset length greater than size in ",
997  "function vector_smallest_index().",exc_einval);
998  }
999  if (k==0 || n==0) {
1000  O2SCL_ERR2("Vector size zero or k zero in ",
1001  "function vector_smallest_index().",exc_einval);
1002  }
1003 
1004  index.resize(k);
1005 
1006  // [GSL] Take the first element
1007  size_t j=1;
1008  data_t xbound=data[0];
1009  index[0]=0;
1010 
1011  // [GSL] Examine the remaining elements
1012  for(size_t i=1;i<n;i++) {
1013  data_t xi=data[i];
1014  if (j<k) {
1015  j++;
1016  } else if (xi>=xbound) {
1017  continue;
1018  }
1019  size_t i1;
1020  for(i1=j-1;i1>0;i1--) {
1021  if (xi>data[index[i1-1]]) break;
1022  index[i1]=index[i1-1];
1023  }
1024  index[i1]=i;
1025  xbound=data[index[j-1]];
1026  }
1027  return;
1028  }
1029 
1030  /** \brief Find the indexes of the k smallest entries of a vector
1031  */
1032  template<class vec_t, class data_t, class vec_size_t>
1033  void vector_smallest_index(const vec_t &data, size_t k,
1034  vec_size_t &index) {
1035  size_t n=data.size();
1036  if (index.size()<k) index.resize(k);
1037  return o2scl::vector_smallest_index<vec_t,data_t,vec_size_t>
1038  (n,data,k,index);
1039  }
1040 
1041  /** \brief Find the k largest entries of the first \c n elements
1042  of a vector
1043 
1044  Given a vector \c data of size \c n this sets the first \c k
1045  entries of the vector \c largest to the k largest entries from
1046  vector \c data in descending order. The vector \c largest must
1047  be allocated beforehand to hold at least \c k elements.
1048 
1049  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1050 
1051  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1052  \f$ k << N \f$.
1053 
1054  If \c k is zero, then this function does nothing and
1055  returns \ref o2scl::success .
1056  */
1057  template<class vec_t, class data_t>
1058  void vector_largest(size_t n, vec_t &data, size_t k, vec_t &largest) {
1059  if (k>n) {
1060  O2SCL_ERR2("Subset length greater than size in ",
1061  "function vector_largest().",exc_einval);
1062  }
1063  if (k==0 || n==0) {
1064  O2SCL_ERR2("Vector size zero or k zero in ",
1065  "function vector_largest().",exc_einval);
1066  }
1067 
1068  // Take the first element
1069  size_t j=1;
1070  data_t xbound=data[0];
1071  largest[0]=xbound;
1072 
1073  // Examine the remaining elements
1074  for(size_t i=1;i<n;i++) {
1075  data_t xi=data[i];
1076  if (j<k) {
1077  j++;
1078  } else if (xi<=xbound) {
1079  continue;
1080  }
1081  size_t i1;
1082  for(i1=j-1;i1>0;i1--) {
1083  if (xi<largest[i1-1]) break;
1084  largest[i1]=largest[i1-1];
1085  }
1086  largest[i1]=xi;
1087  xbound=largest[j-1];
1088  }
1089  return;
1090  }
1091 
1092  /** \brief Find the k largest entries of a vector
1093  of a vector
1094 
1095  Given a vector \c data, this function sets the first \c k
1096  entries of the vector \c largest to the k largest entries from
1097  vector \c data in ascending order. The vector \c largest
1098  is resized if necessary to hold at least \c k elements.
1099 
1100  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1101 
1102  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1103  \f$ k << N \f$.
1104 
1105  If \c k is zero, then this function does nothing and
1106  returns \ref o2scl::success .
1107  */
1108  template<class vec_t, class data_t>
1109  void vector_largest(vec_t &data, size_t k, vec_t &largest) {
1110  size_t n=data.size();
1111  if (largest.size()<k) largest.resize(k);
1112  return vector_largest(n,data,k,largest);
1113  }
1114  //@}
1115 
1116  /// \name Vector minimum and maximum functions in src/base/vector.h
1117  //@{
1118  /** \brief Compute the maximum of the first \c n elements of a vector
1119  */
1120  template<class vec_t, class data_t>
1121  data_t vector_max_value(size_t n, const vec_t &data) {
1122 
1123  if (n==0) {
1124  O2SCL_ERR("Sent size=0 to vector_max_value().",exc_efailed);
1125  }
1126  data_t max=data[0];
1127  for(size_t i=1;i<n;i++) {
1128  if (data[i]>max) {
1129  max=data[i];
1130  }
1131  }
1132  return max;
1133  }
1134 
1135  /** \brief Compute the maximum value of a vector
1136  */
1137  template<class vec_t, class data_t>
1138  data_t vector_max_value(const vec_t &data) {
1139 
1140  size_t n=data.size();
1141  if (n==0) {
1142  O2SCL_ERR("Sent empty vector to vector_max_value().",exc_efailed);
1143  }
1144  data_t max=data[0];
1145  for(size_t i=1;i<n;i++) {
1146  if (data[i]>max) {
1147  max=data[i];
1148  }
1149  }
1150  return max;
1151  }
1152 
1153  /** \brief Compute the index which holds the
1154  maximum of the first \c n elements of a vector
1155  */
1156  template<class vec_t, class data_t>
1157  size_t vector_max_index(size_t n, const vec_t &data) {
1158 
1159  if (n==0) {
1160  O2SCL_ERR("Sent size=0 to vector_max_index().",exc_efailed);
1161  }
1162  data_t max=data[0];
1163  size_t ix=0;
1164  for(size_t i=1;i<n;i++) {
1165  if (data[i]>max) {
1166  max=data[i];
1167  ix=i;
1168  }
1169  }
1170  return ix;
1171  }
1172 
1173  /** \brief Compute the maximum of the first \c n elements of a vector
1174  */
1175  template<class vec_t, class data_t>
1176  void vector_max(size_t n, const vec_t &data, size_t &index,
1177  data_t &val) {
1178 
1179  if (n==0) {
1180  O2SCL_ERR("Sent size=0 to vector_max().",exc_efailed);
1181  }
1182  val=data[0];
1183  index=0;
1184  for(size_t i=1;i<n;i++) {
1185  if (data[i]>val) {
1186  val=data[i];
1187  index=i;
1188  }
1189  }
1190  return;
1191  }
1192 
1193  /** \brief Compute the minimum of the first \c n elements of a vector
1194  */
1195  template<class vec_t, class data_t>
1196  data_t vector_min_value(size_t n, const vec_t &data) {
1197 
1198  if (n==0) {
1199  O2SCL_ERR("Sent size=0 to vector_min_value().",exc_efailed);
1200  }
1201  data_t min=data[0];
1202  for(size_t i=1;i<n;i++) {
1203  if (data[i]<min) {
1204  min=data[i];
1205  }
1206  }
1207  return min;
1208  }
1209 
1210  /** \brief Compute the minimum value in a vector
1211  */
1212  template<class vec_t, class data_t>
1213  data_t vector_min_value(const vec_t &data) {
1214 
1215  size_t n=data.size();
1216  if (n==0) {
1217  O2SCL_ERR("Sent empty vector to vector_min_value().",exc_efailed);
1218  }
1219  data_t min=data[0];
1220  for(size_t i=1;i<n;i++) {
1221  if (data[i]<min) {
1222  min=data[i];
1223  }
1224  }
1225  return min;
1226  }
1227 
1228  /** \brief Compute the index which holds the
1229  minimum of the first \c n elements of a vector
1230  */
1231  template<class vec_t, class data_t>
1232  size_t vector_min_index(size_t n, const vec_t &data) {
1233 
1234  if (n==0) {
1235  O2SCL_ERR("Sent size=0 to vector_min_index().",exc_efailed);
1236  }
1237  data_t min=data[0];
1238  size_t ix=0;
1239  for(size_t i=1;i<n;i++) {
1240  if (data[i]<min) {
1241  min=data[i];
1242  ix=i;
1243  }
1244  }
1245  return ix;
1246  }
1247 
1248  /** \brief Compute the minimum of the first \c n elements of a vector
1249  */
1250  template<class vec_t, class data_t>
1251  void vector_min(size_t n, const vec_t &data, size_t &index,
1252  data_t &val) {
1253 
1254  if (n==0) {
1255  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1256  }
1257  val=data[0];
1258  index=0;
1259  for(size_t i=1;i<n;i++) {
1260  if (data[i]<val) {
1261  val=data[i];
1262  index=i;
1263  }
1264  }
1265  return;
1266  }
1267 
1268  /** \brief Compute the minimum and maximum of the first
1269  \c n elements of a vector
1270  */
1271  template<class vec_t, class data_t>
1272  void vector_minmax_value(size_t n, vec_t &data,
1273  data_t &min, data_t &max) {
1274 
1275  if (n==0) {
1276  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1277  }
1278  min=data[0];
1279  max=min;
1280  for(size_t i=1;i<n;i++) {
1281  if (data[i]<min) {
1282  min=data[i];
1283  }
1284  if (data[i]>max) {
1285  max=data[i];
1286  }
1287  }
1288  return;
1289  }
1290 
1291  /** \brief Compute the minimum and maximum of the first
1292  \c n elements of a vector
1293  */
1294  template<class vec_t, class data_t>
1295  void vector_minmax_index(size_t n, vec_t &data,
1296  size_t &ix_min, size_t &ix_max) {
1297 
1298  if (n==0) {
1299  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1300  }
1301  data_t min=data[0];
1302  data_t max=min;
1303  ix_min=0;
1304  ix_max=0;
1305  for(size_t i=1;i<n;i++) {
1306  if (data[i]<min) {
1307  min=data[i];
1308  ix_min=i;
1309  }
1310  if (data[i]>max) {
1311  max=data[i];
1312  ix_max=i;
1313  }
1314  }
1315  return;
1316  }
1317 
1318  /** \brief Compute the minimum and maximum of the first
1319  \c n elements of a vector
1320  */
1321  template<class vec_t, class data_t>
1322  void vector_minmax(size_t n, vec_t &data,
1323  size_t &ix_min, data_t &min,
1324  size_t &ix_max, data_t &max) {
1325 
1326  if (n==0) {
1327  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1328  }
1329  min=data[0];
1330  max=min;
1331  ix_min=0;
1332  ix_max=0;
1333  for(size_t i=1;i<n;i++) {
1334  if (data[i]<min) {
1335  min=data[i];
1336  ix_min=i;
1337  }
1338  if (data[i]>max) {
1339  max=data[i];
1340  ix_max=i;
1341  }
1342  }
1343  return;
1344  }
1345  //@}
1346 
1347  /// \name Extrema of vectors through quadratic fit in src/base/vector.h
1348  //@{
1349  /** \brief Maximum of vector by quadratic fit
1350  */
1351  template<class vec_t, class data_t>
1352  data_t vector_max_quad(size_t n, const vec_t &data) {
1353  size_t ix=vector_max_index<vec_t,data_t>(n,data);
1354  if (ix==0) {
1355  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1356  } else if (ix==n-1) {
1357  return quadratic_extremum_y<data_t>
1358  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1359  }
1360  return quadratic_extremum_y<data_t>
1361  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1362  }
1363 
1364  /** \brief Maximum of vector by quadratic fit
1365  */
1366  template<class vec_t, class data_t>
1367  data_t vector_max_quad(size_t n, const vec_t &x, const vec_t &y) {
1368  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1369  if (ix==0 || ix==n-1) return y[ix];
1370  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1371  y[ix-1],y[ix],y[ix+1]);
1372  }
1373 
1374  /** \brief Location of vector maximum by quadratic fit
1375  */
1376  template<class vec_t, class data_t>
1377  data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1378  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1379  if (ix==0 || ix==n-1) return y[ix];
1380  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1381  y[ix-1],y[ix],y[ix+1]);
1382  }
1383 
1384  /** \brief Minimum of vector by quadratic fit
1385  */
1386  template<class vec_t, class data_t>
1387  data_t vector_min_quad(size_t n, const vec_t &data) {
1388  size_t ix=vector_min_index<vec_t,data_t>(n,data);
1389  if (ix==0) {
1390  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1391  } else if (ix==n-1) {
1392  return quadratic_extremum_y<data_t>
1393  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1394  }
1395  return quadratic_extremum_y<data_t>
1396  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1397  }
1398 
1399  /** \brief Minimum of vector by quadratic fit
1400  */
1401  template<class vec_t, class data_t>
1402  data_t vector_min_quad(size_t n, const vec_t &x, const vec_t &y) {
1403  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1404  if (ix==0 || ix==n-1) return y[ix];
1405  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1406  y[ix-1],y[ix],y[ix+1]);
1407  }
1408 
1409  /** \brief Location of vector minimum by quadratic fit
1410  */
1411  template<class vec_t, class data_t>
1412  data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1413  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1414  if (ix==0 || ix==n-1) return y[ix];
1415  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1416  y[ix-1],y[ix],y[ix+1]);
1417  }
1418  //@}
1419 
1420  /// \name Matrix minimum and maximum functions in src/base/vector.h
1421  //@{
1422  /** \brief Compute the maximum of the lower-left part of a matrix
1423  */
1424  template<class mat_t, class data_t>
1425  data_t matrix_max_value(size_t m, const size_t n, const mat_t &data) {
1426 
1427  if (m==0 || n==0) {
1428  std::string str=((std::string)"Matrix with zero size (")+
1429  o2scl::itos(m)+","+o2scl::itos(n)+") in "+
1430  "matrix_max_value().";
1431  O2SCL_ERR(str.c_str(),exc_einval);
1432  }
1433  data_t max=data(0,0);
1434  for(size_t i=0;i<m;i++) {
1435  for(size_t j=0;j<n;j++) {
1436  if (data(i,j)>max) {
1437  max=data(i,j);
1438  }
1439  }
1440  }
1441  return max;
1442  }
1443 
1444  /** \brief Compute the maximum of a matrix
1445  */
1446  template<class mat_t, class data_t> data_t
1447  matrix_max_value(const mat_t &data) {
1448  size_t m=data.size1();
1449  size_t n=data.size2();
1450  if (m==0 || n==0) {
1451  std::string str=((std::string)"Matrix with zero size (")+
1452  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1453  "matrix_max_value().";
1454  O2SCL_ERR(str.c_str(),exc_einval);
1455  }
1456  data_t max=data(0,0);
1457  for(size_t i=0;i<m;i++) {
1458  for(size_t j=0;j<n;j++) {
1459  if (data(i,j)>max) {
1460  max=data(i,j);
1461  }
1462  }
1463  }
1464  return max;
1465  }
1466 
1467  /** \brief Compute the maximum of a matrix
1468  */
1469  template<class mat_t> double
1470  matrix_max_value_double(const mat_t &data) {
1471  size_t m=data.size1();
1472  size_t n=data.size2();
1473  if (m==0 || n==0) {
1474  std::string str=((std::string)"Matrix with zero size (")+
1475  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1476  "matrix_max_value_double().";
1477  O2SCL_ERR(str.c_str(),exc_einval);
1478  }
1479  double max=data(0,0);
1480  for(size_t i=0;i<m;i++) {
1481  for(size_t j=0;j<n;j++) {
1482  if (data(i,j)>max) {
1483  max=data(i,j);
1484  }
1485  }
1486  }
1487  return max;
1488  }
1489 
1490  /** \brief Compute the maximum of a matrix and return
1491  the indices of the maximum element
1492  */
1493  template<class mat_t, class data_t>
1494  void matrix_max_index(size_t m, size_t n, const mat_t &data,
1495  size_t &i_max, size_t &j_max, data_t &max) {
1496 
1497  if (m==0 || n==0) {
1498  std::string str=((std::string)"Matrix with zero size (")+
1499  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1500  "matrix_max_index().";
1501  O2SCL_ERR(str.c_str(),exc_einval);
1502  }
1503  max=data(0,0);
1504  i_max=0;
1505  j_max=0;
1506  for(size_t i=0;i<m;i++) {
1507  for(size_t j=0;j<n;j++) {
1508  if (data(i,j)>max) {
1509  max=data(i,j);
1510  i_max=i;
1511  j_max=j;
1512  }
1513  }
1514  }
1515  return;
1516  }
1517 
1518  /** \brief Compute the maximum of a matrix and return
1519  the indices of the maximum element
1520  */
1521  template<class mat_t, class data_t>
1522  void matrix_max_index(const mat_t &data,
1523  size_t &i_max, size_t &j_max, data_t &max) {
1524 
1525  size_t m=data.size1();
1526  size_t n=data.size2();
1527  if (m==0 || n==0) {
1528  std::string str=((std::string)"Matrix with zero size (")+
1529  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1530  "matrix_max_index().";
1531  O2SCL_ERR(str.c_str(),exc_einval);
1532  }
1533  max=data(0,0);
1534  i_max=0;
1535  j_max=0;
1536  for(size_t i=0;i<m;i++) {
1537  for(size_t j=0;j<n;j++) {
1538  if (data(i,j)>max) {
1539  max=data(i,j);
1540  i_max=i;
1541  j_max=j;
1542  }
1543  }
1544  }
1545  return;
1546  }
1547 
1548  /** \brief Compute the minimum of a matrix
1549  */
1550  template<class mat_t, class data_t>
1551  data_t matrix_min_value(size_t m, size_t n, const mat_t &data) {
1552 
1553  if (m==0 || n==0) {
1554  std::string str=((std::string)"Matrix with zero size (")+
1555  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1556  "matrix_min_value().";
1557  O2SCL_ERR(str.c_str(),exc_einval);
1558  }
1559  data_t min=data(0,0);
1560  for(size_t i=0;i<m;i++) {
1561  for(size_t j=0;j<n;j++) {
1562  if (data(i,j)<min) {
1563  min=data(i,j);
1564  }
1565  }
1566  }
1567  return min;
1568  }
1569 
1570  /** \brief Compute the minimum of a matrix
1571  */
1572  template<class mat_t, class data_t>
1573  data_t matrix_min_value(const mat_t &data) {
1574 
1575  size_t m=data.size1();
1576  size_t n=data.size2();
1577  if (m==0 || n==0) {
1578  std::string str=((std::string)"Matrix with zero size (")+
1579  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1580  "matrix_min_value().";
1581  O2SCL_ERR(str.c_str(),exc_einval);
1582  }
1583  data_t min=data(0,0);
1584  for(size_t i=0;i<m;i++) {
1585  for(size_t j=0;j<n;j++) {
1586  if (data(i,j)<min) {
1587  min=data(i,j);
1588  }
1589  }
1590  }
1591  return min;
1592  }
1593 
1594  /** \brief Compute the minimum of a matrix
1595  */
1596  template<class mat_t>
1597  double matrix_min_value_double(const mat_t &data) {
1598 
1599  size_t m=data.size1();
1600  size_t n=data.size2();
1601  if (m==0 || n==0) {
1602  std::string str=((std::string)"Matrix with zero size (")+
1603  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1604  "matrix_min_value().";
1605  O2SCL_ERR(str.c_str(),exc_einval);
1606  }
1607  double min=data(0,0);
1608  for(size_t i=0;i<m;i++) {
1609  for(size_t j=0;j<n;j++) {
1610  if (data(i,j)<min) {
1611  min=data(i,j);
1612  }
1613  }
1614  }
1615  return min;
1616  }
1617 
1618  /** \brief Compute the minimum of a matrix and return
1619  the indices of the minimum element
1620  */
1621  template<class mat_t, class data_t>
1622  void matrix_min_index(size_t n, size_t m, const mat_t &data,
1623  size_t &i_min, size_t &j_min, data_t &min) {
1624 
1625  if (m==0 || n==0) {
1626  std::string str=((std::string)"Matrix with zero size (")+
1627  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1628  "matrix_min_index().";
1629  O2SCL_ERR(str.c_str(),exc_einval);
1630  }
1631  min=data(0,0);
1632  i_min=0;
1633  j_min=0;
1634  for(size_t i=0;i<m;i++) {
1635  for(size_t j=0;j<n;j++) {
1636  if (data(i,j)<min) {
1637  min=data(i,j);
1638  i_min=i;
1639  j_min=j;
1640  }
1641  }
1642  }
1643  return;
1644  }
1645 
1646  /** \brief Compute the minimum of a matrix and return
1647  the indices of the minimum element
1648  */
1649  template<class mat_t, class data_t>
1650  void matrix_min_index(const mat_t &data,
1651  size_t &i_min, size_t &j_min, data_t &min) {
1652 
1653  size_t m=data.size1();
1654  size_t n=data.size2();
1655  if (m==0 || n==0) {
1656  std::string str=((std::string)"Matrix with zero size (")+
1657  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1658  "matrix_min_index().";
1659  O2SCL_ERR(str.c_str(),exc_einval);
1660  }
1661  min=data(0,0);
1662  i_min=0;
1663  j_min=0;
1664  for(size_t i=0;i<m;i++) {
1665  for(size_t j=0;j<n;j++) {
1666  if (data(i,j)<min) {
1667  min=data(i,j);
1668  i_min=i;
1669  j_min=j;
1670  }
1671  }
1672  }
1673  return;
1674  }
1675 
1676  /** \brief Compute the minimum and maximum of a matrix
1677  */
1678  template<class mat_t, class data_t>
1679  void matrix_minmax(size_t n, size_t m, const mat_t &data,
1680  data_t &min, data_t &max) {
1681 
1682  if (n==0 || m==0) {
1683  O2SCL_ERR("Sent size=0 to matrix_min().",exc_efailed);
1684  }
1685  min=data(0,0);
1686  max=data(0,0);
1687  for(size_t i=0;i<n;i++) {
1688  for(size_t j=0;j<m;j++) {
1689  if (data(i,j)<min) {
1690  min=data(i,j);
1691  } else if (data(i,j)>max) {
1692  max=data(i,j);
1693  }
1694  }
1695  }
1696  return;
1697  }
1698 
1699  /** \brief Compute the minimum and maximum of a matrix
1700  */
1701  template<class mat_t, class data_t>
1702  void matrix_minmax(const mat_t &data,
1703  data_t &min, data_t &max) {
1704  return matrix_minmax<mat_t,data_t>
1705  (data.size1(),data.size2(),data,min,max);
1706  }
1707 
1708  /** \brief Compute the minimum and maximum of a matrix and
1709  return their locations
1710  */
1711  template<class mat_t, class data_t>
1712  void matrix_minmax_index(size_t n, size_t m, const mat_t &data,
1713  size_t &i_min, size_t &j_min, data_t &min,
1714  size_t &i_max, size_t &j_max, data_t &max) {
1715 
1716  if (n==0 || m==0) {
1717  O2SCL_ERR2("Sent size=0 to function ",
1718  "matrix_minmax_index().",exc_efailed);
1719  }
1720  min=data(0,0);
1721  i_min=0;
1722  j_min=0;
1723  max=data(0,0);
1724  i_max=0;
1725  j_max=0;
1726  for(size_t i=0;i<n;i++) {
1727  for(size_t j=0;j<m;j++) {
1728  if (data(i,j)<min) {
1729  min=data(i,j);
1730  i_min=i;
1731  j_min=j;
1732  } else if (data(i,j)>max) {
1733  max=data(i,j);
1734  i_max=i;
1735  j_max=j;
1736  }
1737  }
1738  }
1739  return;
1740  }
1741 
1742  /** \brief Compute the sum of matrix elements
1743  */
1744  template<class mat_t, class data_t>
1745  data_t matrix_sum(size_t m, size_t n, const mat_t &data) {
1746 
1747  data_t sum=0.0;
1748  for(size_t i=0;i<m;i++) {
1749  for(size_t j=0;j<n;j++) {
1750  sum+=data(i,j);
1751  }
1752  }
1753  return sum;
1754  }
1755 
1756  /** \brief Compute the sum of matrix elements
1757  */
1758  template<class mat_t, class data_t>
1759  data_t matrix_sum(const mat_t &data) {
1760  return matrix_sum<mat_t,data_t>(data.size1(),data.size2(),data);
1761  }
1762  //@}
1763 
1764  /// \name Searching vectors and matrices in src/base/vector.h
1765  //@{
1766  /** \brief Lookup the value \c x0 in the first \c n elements of
1767  vector \c x
1768 
1769  The function finds the element among the first \c n elements of
1770  \c x which is closest to the value \c x0. It ignores all
1771  elements in \c x which are not finite. If the vector is empty,
1772  or if all of the first \c n elements in \c x are not finite,
1773  then the error handler will be called.
1774 
1775  This function works for all classes \c vec_t where an operator[]
1776  is defined which returns a double (either as a value or a
1777  reference).
1778  */
1779  template<class vec_t>
1780  size_t vector_lookup(size_t n, const vec_t &x, double x0) {
1781  if (n==0) {
1782  O2SCL_ERR("Empty vector in function vector_lookup().",
1783  exc_einval);
1784  return 0;
1785  }
1786  size_t row=0, i=0;
1787  while(!std::isfinite(x[i]) && i<n-1) i++;
1788  if (i==n-1) {
1789  O2SCL_ERR2("Entire vector not finite in ",
1790  "function vector_lookup()",exc_einval);
1791  return 0;
1792  }
1793  double best=x[i], bdiff=fabs(x[i]-x0);
1794  for(;i<n;i++) {
1795  if (std::isfinite(x[i]) && fabs(x[i]-x0)<bdiff) {
1796  row=i;
1797  best=x[i];
1798  bdiff=fabs(x[i]-x0);
1799  }
1800  }
1801  return row;
1802  }
1803 
1804  /** \brief Lookup element \c x0 in vector \c x
1805 
1806  This function finds the element in vector \c x which is closest
1807  to \c x0. It ignores all elements in \c x which are not finite.
1808  If the vector is empty, or if all of the
1809  elements in \c x are not finite, then the error handler will be
1810  called.
1811 
1812  This function works for all classes \c vec_t with a
1813  <tt>size()</tt> method and where an operator[] is defined which
1814  returns a double (either as a value or a reference).
1815  */
1816  template<class vec_t>
1817  size_t vector_lookup(const vec_t &x, double x0) {
1818  return vector_lookup(x.size(),x,x0);
1819  }
1820 
1821  /** \brief Lookup an element in the first $(m,n)$ entries in a matrix
1822 
1823  Return the location <tt>(i,j)</tt> of the element closest to
1824  \c x0.
1825  */
1826  template<class mat_t>
1827  void matrix_lookup(size_t m, size_t n, const mat_t &A,
1828  double x0, size_t &i, size_t &j) {
1829  if (m==0 || n==0) {
1830  O2SCL_ERR("Empty matrix in matrix_lookup().",
1831  exc_einval);
1832  }
1833  double dist=0.0;
1834  bool found_one=false;
1835  for(size_t i2=0;i2<m;i2++) {
1836  for(size_t j2=0;j2<n;j2++) {
1837  if (std::isfinite(A(i,j))) {
1838  if (found_one==false) {
1839  dist=fabs(A(i,j)-x0);
1840  found_one=true;
1841  i=i2;
1842  j=j2;
1843  } else {
1844  if (fabs(A(i,j)-x0)<dist) {
1845  dist=fabs(A(i,j)-x0);
1846  i=i2;
1847  j=j2;
1848  }
1849  }
1850  }
1851  }
1852  }
1853  if (found_one==false) {
1854  O2SCL_ERR2("Entire matrix not finite in ",
1855  "function matrix_lookup()",exc_einval);
1856  }
1857  return;
1858  }
1859 
1860  /** \brief Lookup an element in a matrix
1861 
1862  Return the location <tt>(i,j)</tt> of the element closest to
1863  \c x0.
1864  */
1865  template<class mat_t>
1866  void matrix_lookup(const mat_t &A,
1867  double x0, size_t &i, size_t &j) {
1868  matrix_lookup(A.size1(),A.size2(),x0,i,j);
1869  return;
1870  }
1871 
1872  /** \brief Binary search a part of an increasing vector for
1873  <tt>x0</tt>.
1874 
1875  This function performs a binary search of between
1876  <tt>x[lo]</tt> and <tt>x[hi]</tt>. It returns
1877  - \c lo if \c x0 < <tt>x[lo]</tt>
1878  - \c i if <tt>x[i]</tt> <= \c x0 < <tt>x[i+2]</tt>
1879  for \c lo <= \c i < \c hi
1880  - \c hi-1 if \c x0 >= \c <tt>x[hi-1]</tt>
1881 
1882  This function is designed to find the interval containing \c x0,
1883  not the index of the element closest to x0. To perform the
1884  latter operation, you can use \ref vector_lookup().
1885 
1886  The element at <tt>x[hi]</tt> is never referenced by this
1887  function. The parameter \c hi can be either the index of the
1888  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1889  with starting index <tt>0</tt>), or the index of one element
1890  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1891  index <tt>0</tt>) for a depending on whether or not the user
1892  wants to allow the function to return the index of the last
1893  element.
1894 
1895  This function operates in the same way as
1896  <tt>gsl_interp_bsearch()</tt>.
1897 
1898  The operation of this function is undefined if the data is
1899  not strictly monotonic, i.e. if some of the data elements are
1900  equal.
1901 
1902  This function will call the error handler if \c lo is
1903  greater than \c hi.
1904  */
1905  template<class vec_t, class data_t>
1906  size_t vector_bsearch_inc(const data_t x0, const vec_t &x,
1907  size_t lo, size_t hi) {
1908  if (lo>hi) {
1909  O2SCL_ERR2("Low and high indexes backwards in ",
1910  "function vector_bsearch_inc().",exc_einval);
1911  }
1912  while (hi>lo+1) {
1913  size_t i=(hi+lo)/2;
1914  if (x[i]>x0) {
1915  hi=i;
1916  } else {
1917  lo=i;
1918  }
1919  }
1920 
1921  return lo;
1922  }
1923 
1924  /** \brief Binary search a part of an decreasing vector for
1925  <tt>x0</tt>.
1926 
1927  This function performs a binary search of between
1928  <tt>x[lo]</tt> and <tt>x[hi]</tt> (inclusive). It returns
1929  - \c lo if \c x0 > <tt>x[lo]</tt>
1930  - \c i if <tt>x[i]</tt> >= \c x0 > <tt>x[i+1]</tt>
1931  for \c lo <= \c i < \c hi
1932  - \c hi-1 if \c x0 <= \c <tt>x[hi-1]</tt>
1933 
1934  The element at <tt>x[hi]</tt> is never referenced by this
1935  function. The parameter \c hi can be either the index of the
1936  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1937  with starting index <tt>0</tt>), or the index of one element
1938  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1939  index <tt>0</tt>) for a depending on whether or not the user
1940  wants to allow the function to return the index of the last
1941  element.
1942 
1943  The operation of this function is undefined if the data is
1944  not strictly monotonic, i.e. if some of the data elements are
1945  equal.
1946 
1947  This function will call the error handler if \c lo is
1948  greater than \c hi.
1949  */
1950  template<class vec_t, class data_t>
1951  size_t vector_bsearch_dec(const data_t x0, const vec_t &x,
1952  size_t lo, size_t hi) {
1953  if (lo>hi) {
1954  O2SCL_ERR2("Low and high indexes backwards in ",
1955  "function vector_bsearch_dec().",exc_einval);
1956  }
1957  while (hi>lo+1) {
1958  size_t i=(hi+lo)/2;
1959  if (x[i]<x0) {
1960  hi=i;
1961  } else {
1962  lo=i;
1963  }
1964  }
1965 
1966  return lo;
1967  }
1968 
1969  /** \brief Binary search a part of a monotonic vector for
1970  <tt>x0</tt>.
1971 
1972  This wrapper just calls \ref o2scl::vector_bsearch_inc() or
1973  \ref o2scl::vector_bsearch_dec() depending on the ordering
1974  of \c x.
1975  */
1976  template<class vec_t, class data_t>
1977  size_t vector_bsearch(const data_t x0, const vec_t &x,
1978  size_t lo, size_t hi) {
1979  if (x[lo]<x[hi-1]) {
1980  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1981  }
1982  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
1983  }
1984 
1985  /** \brief Binary search a monotonic vector for
1986  <tt>x0</tt>.
1987 
1988  This function calls \ref o2scl::vector_bsearch_inc() or
1989  \ref o2scl::vector_bsearch_dec() depending on the ordering
1990  of \c x.
1991  */
1992  template<class vec_t, class data_t>
1993  size_t vector_bsearch(const data_t x0, const vec_t &x) {
1994  size_t lo=0;
1995  size_t hi=x.size();
1996  if (x[lo]<x[hi-1]) {
1997  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1998  }
1999  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
2000  }
2001  //@}
2002 
2003  /// \name Ordering and finite tests in src/base/vector.h
2004  //@{
2005  /** \brief Test if the first \c n elements of a vector are
2006  monotonic and increasing or decreasing
2007 
2008  If \c n is zero or one, this function will return 0 without
2009  calling the error handler. If all the vector's elements are equal,
2010  this function will return 3. Otherwise, if the vector is not
2011  monotonic, then this function will return 0. Finally, if the
2012  vector is nondecreasing (increasing or equal intervals), this
2013  function will return 1, and if the vector is nonincreasing
2014  (decreasing or equal intervals), this function will return 2.
2015  This function assumes that simple comparison operators have been
2016  defined for the type of each vector element.
2017  */
2018  template<class vec_t>
2019  int vector_is_monotonic(size_t n, vec_t &data) {
2020 
2021  if (n<1) return 0;
2022  if (n<2) {
2023  if (data[0]==data[1]) {
2024  return 3;
2025  } else if (data[0]<data[1]) {
2026  return 1;
2027  } else {
2028  return 2;
2029  }
2030  }
2031 
2032  // Find first non-flat interval
2033  size_t start=0;
2034  bool done=false;
2035  for(size_t i=0;i<n-1 && done==false;i++) {
2036  if (data[i]!=data[i+1]) {
2037  done=true;
2038  } else {
2039  start++;
2040  }
2041  }
2042 
2043  // If all elements in the vector are equal, or if the only
2044  // distinct element is at the end, then return true.
2045  if (done==false) {
2046  return 3;
2047  }
2048 
2049  if (start==n-2) {
2050  if (data[start]<data[start+1]) return 1;
2051  else return 2;
2052  }
2053 
2054  // Determine if the vector is increasing (inc=true) or decreasing
2055  // (inc=false)
2056  bool inc=true;
2057  if (data[start]>data[start+1]) inc=false;
2058 
2059  if (inc) {
2060  for(size_t i=start+1;i<n-1;i++) {
2061  // If there is one decreasing interval, return false
2062  if (data[i]>data[i+1]) return 0;
2063  }
2064  return 1;
2065  }
2066 
2067  // If there is one increasing interval, return false
2068  for(size_t i=start+1;i<n-1;i++) {
2069  if (data[i]<data[i+1]) return 0;
2070  }
2071  return 2;
2072  }
2073 
2074  /** \brief Test if the first \c n elements of a vector are
2075  monotonic and increasing or decreasing
2076 
2077  If \c n is zero or one, this function will return 0 without
2078  calling the error handler. If all the vector's elements are equal,
2079  this function will return 3. Otherwise, if the vector is not
2080  monotonic, then this function will return 0. Finally, if the
2081  vector is nondecreasing (increasing or equal intervals), this
2082  function will return 1, and if the vector is nonincreasing
2083  (decreasing or equal intervals), this function will return 2.
2084  This function assumes that simple comparison operators have been
2085  defined for the type of each vector element.
2086  */
2087  template<class vec_t> int vector_is_monotonic(vec_t &data) {
2088  return vector_is_monotonic(data.size(),data);
2089  }
2090 
2091  /** \brief Test if the first \c n elements of a vector are
2092  strictly monotonic and determine if they are increasing or decreasing
2093 
2094  If \c n is zero this function will return 0 without calling the
2095  error handler. Also, if the vector is not monotonic, this
2096  function will return 0. If the vector is strictly
2097  monotonic, then this function will return 1 if it is
2098  increasing and 2 if it is decreasing.
2099  */
2100  template<class vec_t>
2101  int vector_is_strictly_monotonic(size_t n, vec_t &data) {
2102 
2103  if (n<1) return 0;
2104 
2105  // Determine if the vector is increasing (inc=true) or decreasing
2106  // (inc=false)
2107  bool inc=true;
2108  if (data[0]==data[1]) {
2109  return 0;
2110  } else if (data[0]>data[1]) {
2111  inc=false;
2112  }
2113 
2114  if (inc) {
2115  for(size_t i=1;i<n-1;i++) {
2116  // If there is one nonincreasing interval, return 0
2117  if (data[i]>=data[i+1]) return 0;
2118  }
2119  return 1;
2120  }
2121 
2122  // If there is one increasing interval, return 0
2123  for(size_t i=1;i<n-1;i++) {
2124  if (data[i]<=data[i+1]) return 0;
2125  }
2126  return 2;
2127  }
2128 
2129  /** \brief Test if the first \c n elements of a vector are
2130  strictly monotonic and determine if they are increasing or decreasing
2131 
2132  If \c n is zero this function will return 0 without calling the
2133  error handler. Also, if the vector is not monotonic, this
2134  function will return 0. If the vector is strictly
2135  monotonic, then this function will return 1 if it is
2136  increasing and 2 if it is decreasing.
2137  */
2138  template<class vec_t>
2139  int vector_is_strictly_monotonic(vec_t &data) {
2140  return vector_is_strictly_monotonic(data.size(),data);
2141  }
2142 
2143  /** \brief Test if the first \c n elements of a vector are finite
2144 
2145  If \c n is zero, this will return true without throwing
2146  an exception.
2147 
2148  The corresponding tests for matrix functions are
2149  in clbas_base.h .
2150  */
2151  template<class vec_t>
2152  bool vector_is_finite(size_t n, vec_t &data) {
2153  for(size_t i=0;i<n;i++) {
2154  if (!std::isfinite(data[i])) return false;
2155  }
2156  return true;
2157  }
2158 
2159  /** \brief Test if a vector is finite
2160 
2161  If \c n is zero, this will return true without throwing
2162  an exception.
2163 
2164  The corresponding tests for matrix functions are
2165  in clbas_base.h .
2166  */
2167  template<class vec_t> bool vector_is_finite(vec_t &data) {
2168  return vector_is_finite(data.size(),data);
2169  }
2170  //@}
2171 
2172  /// \name Miscellaneous mathematical functions in src/base/vector.h
2173  //@{
2174  /** \brief Compute the sum of the first \c n elements of a vector
2175 
2176  If \c n is zero, this will return 0 without throwing
2177  an exception.
2178  */
2179  template<class vec_t, class data_t>
2180  data_t vector_sum(size_t n, vec_t &data) {
2181 
2182  data_t sum=0.0;
2183  for(size_t i=0;i<n;i++) {
2184  sum+=data[i];
2185  }
2186  return sum;
2187  }
2188 
2189  /** \brief Create a new vector containing the differences between
2190  adjacent entries
2191  */
2192  template<class vec_t, class rvec_t>
2193  void vector_diffs(const vec_t &v_data, rvec_t &v_diffs) {
2194 
2195  size_t n=v_data.size();
2196  v_diffs.resize(n-1);
2197  for(size_t i=0;i<n-1;i++) {
2198  v_diffs[i]=v_data[i+1]-v_data[i];
2199  }
2200 
2201  return;
2202  }
2203 
2204  /** \brief Compute the sum of all the elements of a vector
2205 
2206  If the vector has zero size, this will return 0 without
2207  calling the error handler.
2208  */
2209  template<class vec_t, class data_t> data_t vector_sum(vec_t &data) {
2210  data_t sum=0.0;
2211  for(size_t i=0;i<data.size();i++) {
2212  sum+=data[i];
2213  }
2214  return sum;
2215  }
2216 
2217  /** \brief Compute the sum of the first \c n elements of a vector
2218  of double-precision numbers
2219 
2220  If \c n is zero, this will return 0 without throwing
2221  an exception.
2222  */
2223  template<class vec_t>double vector_sum_double(size_t n, vec_t &data) {
2224  double sum=0.0;
2225  for(size_t i=0;i<n;i++) {
2226  sum+=data[i];
2227  }
2228  return sum;
2229  }
2230 
2231  /** \brief Compute the sum of all the elements of a vector
2232  of double-precision numbers
2233 
2234  If the vector has zero size, this will return 0 without
2235  calling the error handler.
2236  */
2237  template<class vec_t> double vector_sum_double(vec_t &data) {
2238  double sum=0.0;
2239  for(size_t i=0;i<data.size();i++) {
2240  sum+=data[i];
2241  }
2242  return sum;
2243  }
2244 
2245  /** \brief Compute the norm of the first \c n entries of a
2246  vector of floating-point (single or double precision) numbers
2247 
2248  This function is a more generic version of
2249  \ref o2scl_cblas::dnrm2 .
2250  */
2251  template<class vec_t, class data_t>
2252  data_t vector_norm(size_t n, const vec_t &x) {
2253 
2254  data_t scale = 0.0;
2255  data_t ssq = 1.0;
2256 
2257  if (n <= 0) {
2258  return 0.0;
2259  } else if (n == 1) {
2260  return fabs(x[0]);
2261  }
2262 
2263  for (size_t i = 0; i < n; i++) {
2264  const data_t xx = x[i];
2265 
2266  if (xx != 0.0) {
2267  const data_t ax = fabs(xx);
2268 
2269  if (scale < ax) {
2270  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2271  scale = ax;
2272  } else {
2273  ssq += (ax / scale) * (ax / scale);
2274  }
2275  }
2276 
2277  }
2278 
2279  return scale * sqrt(ssq);
2280  }
2281 
2282  /** \brief Compute the norm of a vector of floating-point
2283  (single or double precision) numbers
2284  */
2285  template<class vec_t, class data_t> data_t vector_norm(const vec_t &x) {
2286  return vector_norm<vec_t,data_t>(x.size(),x);
2287  }
2288 
2289  /** \brief Compute the norm of the first \c n entries of a
2290  vector of double precision numbers
2291 
2292  This function is a more generic version of
2293  \ref o2scl_cblas::dnrm2 .
2294  */
2295  template<class vec_t>
2296  double vector_norm_double(size_t n, const vec_t &x) {
2297 
2298  double scale = 0.0;
2299  double ssq = 1.0;
2300 
2301  if (n <= 0) {
2302  return 0.0;
2303  } else if (n == 1) {
2304  return fabs(x[0]);
2305  }
2306 
2307  for (size_t i = 0; i < n; i++) {
2308  const double xx = x[i];
2309 
2310  if (xx != 0.0) {
2311  const double ax = fabs(xx);
2312 
2313  if (scale < ax) {
2314  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2315  scale = ax;
2316  } else {
2317  ssq += (ax / scale) * (ax / scale);
2318  }
2319  }
2320 
2321  }
2322 
2323  return scale * sqrt(ssq);
2324  }
2325 
2326  /** \brief Compute the norm of a vector of double precision numbers
2327  */
2328  template<class vec_t> double vector_norm_double(const vec_t &x) {
2329  return vector_norm_double<vec_t>(x.size(),x);
2330  }
2331  //@}
2332 
2333  /// \name Other vector and matrix functions in src/base/vector.h
2334  //@{
2335  /** \brief Set the first N entries in a vector to a particular value
2336  */
2337  template<class vec_t, class data_t>
2338  void vector_set_all(size_t N, vec_t &src, data_t val) {
2339  for(size_t i=0;i<N;i++) {
2340  src[i]=val;
2341  }
2342  return;
2343  }
2344 
2345  /** \brief Set all entries in a vector to a particular value
2346  */
2347  template<class vec_t, class data_t>
2348  void vector_set_all(vec_t &src, data_t val) {
2349  o2scl::vector_set_all(src.size(),src,val);
2350  return;
2351  }
2352 
2353  /** \brief Set the first (M,N) entries in a matrix to a particular value
2354  */
2355  template<class mat_t, class data_t>
2356  void matrix_set_all(size_t M, size_t N, mat_t &src, data_t val) {
2357  for(size_t i=0;i<M;i++) {
2358  for(size_t j=0;j<N;j++) {
2359  src(i,j)=val;
2360  }
2361  }
2362  return;
2363  }
2364 
2365  /** \brief Set all entries in a matrix to a particular value
2366  */
2367  template<class mat_t, class data_t>
2368  void matrix_set_all(mat_t &src, data_t val) {
2369  o2scl::matrix_set_all(src.size1(),src.size2(),src,val);
2370  return;
2371  }
2372 
2373  /** \brief From a given vector, create a new vector by removing a
2374  specified element
2375 
2376  This funciton is used in \ref o2scl::interp_krige_optim::qual_fun() .
2377  */
2378  template<class vec_t, class vec2_t>
2379  void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest) {
2380  if (src.size()==0) {
2381  O2SCL_ERR("Empty source vector.",o2scl::exc_einval);
2382  }
2383  if (iout>=src.size()) {
2384  O2SCL_ERR("Requested element beyond end.",o2scl::exc_einval);
2385  }
2386  dest.resize(src.size()-1);
2387  size_t j=0;
2388  for(size_t i=0;i<src.size();i++) {
2389  if (i!=iout) {
2390  dest[j]=src[i];
2391  j++;
2392  }
2393  }
2394  return;
2395  }
2396 
2397  /** \brief From a given vector, create a new vector by removing a
2398  specified element
2399 
2400  This funciton is used in \ref o2scl::interp_krige_optim::qual_fun() .
2401  */
2402  template<class vec_t, class vec2_t>
2403  void vector_copy_jackknife(size_t sz, const vec_t &src,
2404  size_t iout, vec2_t &dest) {
2405 
2406  if (sz==0) {
2407  O2SCL_ERR("Empty source vector.",o2scl::exc_einval);
2408  }
2409  if (iout>=sz) {
2410  O2SCL_ERR("Requested element beyond end.",o2scl::exc_einval);
2411  }
2412  dest.resize(sz-1);
2413  size_t j=0;
2414  for(size_t i=0;i<sz;i++) {
2415  if (i!=iout) {
2416  dest[j]=src[i];
2417  j++;
2418  }
2419  }
2420  return;
2421  }
2422 
2423  /** \brief "Rotate" a vector so that the kth element is now the beginning
2424 
2425  This is a generic template function which will work for
2426  any types \c data_t and \c vec_t for which
2427  - \c data_t has an <tt>operator=</tt>
2428  - <tt>vec_t::operator[]</tt> returns a reference
2429  to an object of type \c data_t
2430 
2431  This function is used, for example, in \ref o2scl::pinside.
2432 
2433  \note This function is not the same as a Givens rotation,
2434  which is typically referred to in BLAS routines as <tt>drot()</tt>.
2435  */
2436  template<class vec_t, class data_t>
2437  void vector_rotate(size_t n, vec_t &data, size_t k) {
2438 
2439  data_t *tmp=new data_t[n];
2440  for(size_t i=0;i<n;i++) {
2441  tmp[i]=data[(i+k)%n];
2442  }
2443  for(size_t i=0;i<n;i++) {
2444  data[i]=tmp[i];
2445  }
2446  delete[] tmp;
2447 
2448  return;
2449  }
2450 
2451  /** \brief Reverse the first \c n elements of a vector
2452 
2453  If \c n is zero, this function will silently do nothing.
2454  */
2455  template<class vec_t, class data_t>
2456  void vector_reverse(size_t n, vec_t &data) {
2457  data_t tmp;
2458 
2459  for(size_t i=0;i<n/2;i++) {
2460  tmp=data[n-1-i];
2461  data[n-1-i]=data[i];
2462  data[i]=tmp;
2463  }
2464  return;
2465  }
2466 
2467  /** \brief Reverse a vector
2468 
2469  If the <tt>size()</tt> method returns zero, this function will
2470  silently do nothing.
2471  */
2472  template<class vec_t, class data_t>
2473  void vector_reverse(vec_t &data) {
2474  data_t tmp;
2475  size_t n=data.size();
2476 
2477  for(size_t i=0;i<n/2;i++) {
2478  tmp=data[n-1-i];
2479  data[n-1-i]=data[i];
2480  data[i]=tmp;
2481  }
2482  return;
2483  }
2484 
2485  /** \brief Reverse the first n elements in a vector of double
2486  precision numbers
2487 
2488  If \c n is zero, this function will silently do nothing.
2489  */
2490  template<class vec_t>
2491  void vector_reverse_double(size_t n, vec_t &data) {
2492  double tmp;
2493 
2494  for(size_t i=0;i<n/2;i++) {
2495  tmp=data[n-1-i];
2496  data[n-1-i]=data[i];
2497  data[i]=tmp;
2498  }
2499  return;
2500  }
2501 
2502  /** \brief Reverse a vector of double precision numbers
2503 
2504  If the <tt>size()</tt> method returns zero, this function will
2505  silently do nothing.
2506  */
2507  template<class vec_t> void vector_reverse_double(vec_t &data) {
2508  double tmp;
2509  size_t n=data.size();
2510 
2511  for(size_t i=0;i<n/2;i++) {
2512  tmp=data[n-1-i];
2513  data[n-1-i]=data[i];
2514  data[i]=tmp;
2515  }
2516  return;
2517  }
2518 
2519  /** \brief Trivial index vector
2520 
2521  This object just returns the index whenever an object in the
2522  vector is requested, i.e. <tt>operator[](i)</tt> always returns
2523  \c i.
2524  */
2525  template<class data_t> class vector_index_vector {
2526  public:
2527  data_t operator[](size_t &i) const {
2528  return i;
2529  }
2530  };
2531 
2532  /** \brief Index vector with a size method
2533 
2534  This object just returns the index whenever an object in the
2535  vector is requested, i.e. <tt>operator[](i)</tt> always returns
2536  \c i.
2537  */
2538  template<class data_t> class vector_index_vector_size {
2539 
2540  protected:
2541 
2542  /// The vector size
2543  size_t n;
2544 
2545  public:
2546 
2547  /** \brief Create an index vector with size \c n_
2548  */
2550  n=n_;
2551  }
2552 
2553  /** \brief Obtain the element with index \c i
2554  */
2555  data_t operator[](size_t &i) const {
2556  if (i>=n) {
2557  O2SCL_ERR("Out of bounds.",o2scl::exc_einval);
2558  }
2559  return i;
2560  }
2561 
2562  /** \brief Get the size of the vector
2563  */
2564  size_t size() const {
2565  return n;
2566  }
2567 
2568  /** \brief Resize the index vector
2569  */
2570  void resize(size_t n_) {
2571  n=n_;
2572  }
2573 
2574  };
2575 
2576  /** \brief Construct a row of a matrix
2577 
2578  This class template works with combinations of ublas
2579  <tt>matrix</tt> and <tt>matrix_row</tt> objects,
2580  <tt>arma::mat</tt> and <tt>arma::rowvec</tt>, and
2581  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2582 
2583  \note When calling this function with ublas objects, the
2584  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2585  otherwise some compilers will use argument dependent lookup and
2586  get (justifiably) confused with matrix_row in the ublas
2587  namespace.
2588 
2589  \note The template parameters must be explicitly specified
2590  when calling this template function.
2591  */
2592  template<class mat_t, class mat_row_t>
2593  mat_row_t matrix_row(mat_t &M, size_t row) {
2594  return mat_row_t(M,row);
2595  }
2596 
2597  /** \brief Generic object which represents a row of a matrix
2598 
2599  \note This class is experimental.
2600 
2601  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2602  to construct a row of a matrix object of type
2603  \code
2604  std::function<double &(size_t,size_t)>
2605  \endcode
2606  */
2607  template<class mat_t> class matrix_row_gen {
2608 
2609  protected:
2610 
2611  /// A reference to the original matrix
2612  mat_t &m_;
2613 
2614  /// The selected row
2615  size_t row_;
2616 
2617  public:
2618 
2619  /// Create a row object from row \c row of matrix \c m
2620  matrix_row_gen(mat_t &m, size_t row) : m_(m), row_(row) {
2621  }
2622 
2623  /// Return a reference to the ith column of the selected row
2624  double &operator[](size_t i) {
2625  return m_(row_,i);
2626  }
2627 
2628  /// Return a const reference to the ith column of the selected row
2629  const double &operator[](size_t i) const {
2630  return m_(row_,i);
2631  }
2632  };
2633 
2634  /** \brief Matrix row object with a constructor and resize method
2635 
2636  This is used in \ref o2scl::ode_iv_solve_grid .
2637  */
2638  template<class mat_t> class matrix_row_gen_ctor {
2639 
2640  protected:
2641 
2642  /// A pointer to the matrix
2643  mat_t *mp;
2644 
2645  /// The selected row
2646  size_t row_;
2647 
2648  /// A matrix to point to
2649  mat_t mat;
2650 
2651  public:
2652 
2653  /// Create a row object from row \c row of matrix \c m
2654  matrix_row_gen_ctor(mat_t &m, size_t row) : mp(&m), row_(row) {
2655  }
2656 
2657  /// Create a row object from row \c row of matrix \c m
2658  matrix_row_gen_ctor(size_t n_cols=0) {
2659  if (n_cols==0) {
2660  mp=0;
2661  } else {
2662  mat.resize(1,n_cols);
2663  mp=&mat;
2664  row_=0;
2665  }
2666  }
2667 
2668  /** \brief Resize
2669  */
2670  void resize(size_t n_cols=0) {
2671  if (n_cols==0) {
2672  mp=0;
2673  mat.resize(0,0);
2674  } else {
2675  mat.resize(1,n_cols);
2676  mp=&mat;
2677  row_=0;
2678  }
2679  return;
2680  }
2681 
2682  /** \brief Return size
2683  */
2684  size_t size() const {
2685  if (mp==0) {
2686  return 0;
2687  }
2688  return mp->size2();
2689  }
2690 
2691  /// Return a reference to the ith column of the selected row
2692  double &operator[](size_t i) {
2693  if (mp==0) {
2694  O2SCL_ERR("No matrix in matrix_row_gen_ctor::operator[].",
2696  }
2697  return (*mp)(row_,i);
2698  }
2699 
2700  /// Return a const reference to the ith column of the selected row
2701  const double &operator[](size_t i) const {
2702  if (mp==0) {
2703  O2SCL_ERR("No matrix in matrix_row_gen_ctor::operator[].",
2705  }
2706  return (*mp)(row_,i);
2707  }
2708  };
2709 
2710  /** \brief Construct a view of the transpose of a matrix
2711 
2712  \note This class is experimental.
2713  */
2714  template<class mat_t> class matrix_view_transpose {
2715 
2716  protected:
2717 
2718  /// A reference to the original matrix
2719  mat_t &m_;
2720 
2721  public:
2722 
2723  /// Create a row object from row \c row of matrix \c m
2724  matrix_view_transpose(mat_t &m) : m_(m) {
2725  }
2726 
2727  /// Return a reference to the ith column of the selected row
2728  double &operator()(size_t i, size_t j) {
2729  return m_(j,i);
2730  }
2731 
2732  /// Return a const reference to the ith column of the selected row
2733  const double &operator()(size_t i, size_t j) const {
2734  return m_(j,i);
2735  }
2736 
2737  /** \brief Return the number of rows
2738  */
2739  size_t size1() const {
2740  return m_.size2();
2741  }
2742 
2743  /** \brief Return the number of columns
2744  */
2745  size_t size2() const {
2746  return m_.size1();
2747  }
2748 
2749 
2750  };
2751 
2752  /** \brief Construct a view of a matrix omtting a specified row
2753 
2754  \note This class is experimental.
2755  */
2756  template<class mat_t> class matrix_view_omit_row {
2757 
2758  protected:
2759 
2760  /// A reference to the original matrix
2761  mat_t &m_;
2762 
2763  size_t ro;
2764 
2765  public:
2766 
2767  /// Create
2768  matrix_view_omit_row(mat_t &m, size_t row_omit) : m_(m) {
2769  ro=row_omit;
2770  }
2771 
2772  /// Return a reference
2773  double &operator()(size_t i, size_t j) {
2774  if (i>=ro) {
2775  return m_(i+1,j);
2776  }
2777  return m_(i,j);
2778  }
2779 
2780  /// Return a const reference
2781  const double &operator()(size_t i, size_t j) const {
2782  if (i>=ro) {
2783  return m_(i+1,j);
2784  }
2785  return m_(i,j);
2786  }
2787 
2788  /** \brief Return the number of rows
2789  */
2790  size_t size1() const {
2791  return m_.size1()-1;
2792  }
2793 
2794  /** \brief Return the number of columns
2795  */
2796  size_t size2() const {
2797  return m_.size2();
2798  }
2799 
2800 
2801  };
2802 
2803  /** \brief Construct a view of a matrix omitting one specified column
2804 
2805  \note This class is experimental.
2806  */
2807  template<class mat_t> class matrix_view_omit_column {
2808 
2809  protected:
2810 
2811  /// A reference to the original matrix
2812  mat_t &m_;
2813 
2814  size_t co;
2815 
2816  public:
2817 
2818  /// Create
2819  matrix_view_omit_column(mat_t &m, size_t column_omit) : m_(m) {
2820  co=column_omit;
2821  }
2822 
2823  /// Return a reference
2824  double &operator()(size_t i, size_t j) {
2825  if (j>=co) {
2826  return m_(i,j+1);
2827  }
2828  return m_(i,j);
2829  }
2830 
2831  /// Return a const reference
2832  const double &operator()(size_t i, size_t j) const {
2833  if (j>=co) {
2834  return m_(i,j+1);
2835  }
2836  return m_(i,j);
2837  }
2838 
2839  /** \brief Return the number of rows
2840  */
2841  size_t size1() const {
2842  return m_.size1();
2843  }
2844 
2845  /** \brief Return the number of columns
2846  */
2847  size_t size2() const {
2848  return m_.size2()-1;
2849  }
2850 
2851 
2852  };
2853 
2854  /** \brief Generic object which represents a row of a const matrix
2855 
2856  \note This class is experimental.
2857 
2858  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2859  to construct a row of a matrix object of type
2860  \code
2861  std::function<double &(size_t,size_t)>
2862  \endcode
2863  */
2864  template<class mat_t> class const_matrix_row_gen {
2865 
2866  protected:
2867 
2868  /// A reference to the original matrix
2869  const mat_t &m_;
2870 
2871  /// The selected row
2872  size_t row_;
2873 
2874  public:
2875 
2876  /// Create a row object from row \c row of matrix \c m
2877  const_matrix_row_gen(const mat_t &m, size_t row) : m_(m), row_(row) {
2878  }
2879 
2880  /// Return a const reference to the ith column of the selected row
2881  const double &operator[](size_t i) const {
2882  return m_(row_,i);
2883  }
2884  };
2885 
2886  /** \brief A simple matrix view object
2887  */
2889 
2890  public:
2891 
2892  /** \brief Return a reference to the element at row \c row
2893  and column \c col
2894  */
2895  const double &operator()(size_t row, size_t col) const;
2896 
2897  /** \brief Return the number of rows
2898  */
2899  size_t size1() const;
2900 
2901  /** \brief Return the number of columns
2902  */
2903  size_t size2() const;
2904 
2905  };
2906 
2907  /** \brief A simple matrix view object
2908  */
2909  class matrix_view {
2910 
2911  public:
2912 
2913  /** \brief Return a reference to the element at row \c row
2914  and column \c col
2915  */
2916  const double &operator()(size_t row, size_t col) const;
2917 
2918  /** \brief Return a reference to the element at row \c row
2919  and column \c col
2920  */
2921  double &operator()(size_t row, size_t col);
2922 
2923  /** \brief Return the number of rows
2924  */
2925  size_t size1() const;
2926 
2927  /** \brief Return the number of columns
2928  */
2929  size_t size2() const;
2930 
2931  };
2932 
2933  /** \brief View a o2scl::table object as a matrix
2934 
2935  \note This stores a pointer to the table and the user must ensure
2936  that the pointer is valid with the matrix view is accessed.
2937 
2938  \future It would be nice to store a reference rather than a
2939  pointer, but this causes problems with \ref o2scl::interpm_idw .
2940  */
2941  template<class vec1_t, class vec2_t=std::vector<vec1_t> >
2943 
2944  protected:
2945 
2946  /// Pointer to the table
2947  vec2_t *vvp;
2948 
2949  public:
2950 
2951  /** \brief Swap method
2952  */
2953  friend void swap(matrix_view_vec_vec &t1,
2954  matrix_view_vec_vec &t2) {
2955  /// Just swap the pointer
2956  std::swap(t1.vvp,t2.vvp);
2957  return;
2958  }
2959 
2961  vvp=0;
2962  }
2963 
2964  /** \brief Create a matrix view object from the specified
2965  table and list of rows
2966  */
2967  matrix_view_vec_vec(vec2_t &vv) {
2968  vvp=&vv;
2969  }
2970 
2971  /** \brief Return the number of rows
2972  */
2973  size_t size1() const {
2974  if (vvp==0) return 0;
2975  return vvp->size();
2976  }
2977 
2978  /** \brief Return the number of columns
2979  */
2980  size_t size2() const {
2981  if (vvp==0) return 0;
2982  if (vvp->size()==0) return 0;
2983  return (*vvp)[0].size();
2984  }
2985 
2986  /** \brief Return a reference to the element at row \c row
2987  and column \c col
2988  */
2989  const double &operator()(size_t row, size_t col) const {
2990  if (vvp==0) {
2991  O2SCL_ERR2("Object empty in ",
2992  "matrix_view_vec_vec::operator().",
2994  }
2995  if (row>=vvp->size()) {
2996  O2SCL_ERR2("Row exceeds max in ",
2997  "matrix_view_vec_vec::operator().",
2999  }
3000  if (col>=(*vvp)[row].size()) {
3001  O2SCL_ERR2("Column exceeds max in ",
3002  "matrix_view_vec_vec::operator().",
3004  }
3005  return (*vvp)[row][col];
3006  }
3007 
3008  /** \brief Return a reference to the element at row \c row
3009  and column \c col
3010  */
3011  double &operator()(size_t row, size_t col) {
3012  if (vvp==0) {
3013  O2SCL_ERR2("Object empty in ",
3014  "matrix_view_vec_vec::operator().",
3016  }
3017  if (row>=vvp->size()) {
3018  O2SCL_ERR2("Row exceeds max in ",
3019  "matrix_view_vec_vec::operator().",
3021  }
3022  if (col>=(*vvp)[row].size()) {
3023  std::cout << row << " " << col << " "
3024  << (*vvp)[row].size() << std::endl;
3025  O2SCL_ERR2("Column exceeds max in ",
3026  "matrix_view_vec_vec::operator().",
3028  }
3029  return (*vvp)[row][col];
3030  }
3031 
3032  };
3033 
3034  /** \brief Construct a column of a matrix
3035 
3036  This class template works with combinations of ublas
3037  <tt>matrix</tt> and <tt>matrix_column</tt> objects,
3038  <tt>arma::mat</tt> and <tt>arma::colvec</tt>, and
3039  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
3040 
3041  \note When calling this function with ublas objects, the
3042  namespace prefix <tt>"o2scl::"</tt> must often be specified,
3043  otherwise some compilers will use argument dependent lookup and
3044  get (justifiably) confused with matrix_column in the ublas
3045  namespace.
3046 
3047  \note The template parameters must be explicitly specified
3048  when calling this template function.
3049  */
3050  template<class mat_t, class mat_column_t>
3051  mat_column_t matrix_column(mat_t &M, size_t column) {
3052  return mat_column_t(M,column);
3053  }
3054 
3055  /** \brief Generic object which represents a column of a matrix
3056 
3057  \note This class is experimental.
3058 
3059  The only requirement on the type <tt>mat_t</tt> is that
3060  it must have an operator(size_t,size_t) method which
3061  accesses elements in the matrix.
3062 
3063  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
3064  to construct a row of a matrix object of type
3065  \code
3066  std::function<double &(size_t,size_t)>
3067  \endcode
3068  */
3069  template<class mat_t> class matrix_column_gen {
3070 
3071  protected:
3072 
3073  /// A reference to the original matrix
3074  mat_t &m_;
3075 
3076  /// The selected column
3077  size_t column_;
3078 
3079  public:
3080 
3081  /// Create a column object from column \c column of matrix \c m
3082  matrix_column_gen(mat_t &m, size_t column) : m_(m), column_(column) {
3083  }
3084 
3085  /// Return a reference to the ith row of the selected column
3086  double &operator[](size_t i) {
3087  return m_(i,column_);
3088  }
3089 
3090  /// Return a const reference to the ith row of the selected column
3091  const double &operator[](size_t i) const {
3092  return m_(i,column_);
3093  }
3094 
3095  };
3096 
3097  /** \brief Generic object which represents a column of a const matrix
3098 
3099  \note This class is experimental.
3100 
3101  The only requirement on the type <tt>mat_t</tt> is that
3102  it must have an operator(size_t,size_t) method which
3103  accesses elements in the matrix.
3104 
3105  This class is used in one of
3106  the \ref o2scl::prob_dens_mdim_gaussian constructors.
3107  */
3108  template<class mat_t> class const_matrix_column_gen {
3109 
3110  protected:
3111 
3112  /// A reference to the original matrix
3113  const mat_t &m_;
3114 
3115  /// The selected column
3116  size_t column_;
3117 
3118  public:
3119 
3120  /// Create a column object from column \c column of matrix \c m
3121  const_matrix_column_gen(const mat_t &m, size_t column) :
3122  m_(m), column_(column) {
3123  }
3124 
3125  /// Return a const reference to the ith row of the selected column
3126  const double &operator[](size_t i) const {
3127  return m_(i,column_);
3128  }
3129 
3130  };
3131 
3132  /** \brief Output the first \c n elements of a vector to a stream,
3133  \c os
3134 
3135  No trailing space is output after the last element, and an
3136  endline is output only if \c endline is set to \c true. If the
3137  parameter \c n is zero, this function silently does nothing.
3138 
3139  This works with any class <tt>vec_t</tt> which has an operator[]
3140  which returns either the value of or a reference to the ith
3141  element and the element type has its own output operator which
3142  has been defined.
3143  */
3144  template<class vec_t>
3145  void vector_out(std::ostream &os, size_t n, const vec_t &v,
3146  bool endline=false) {
3147 
3148  // This next line is important since n-1 is not well-defined if n=0
3149  if (n==0) {
3150  if (endline) os << std::endl;
3151  return;
3152  }
3153 
3154  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
3155  os << v[n-1];
3156  if (endline) os << std::endl;
3157  return;
3158  }
3159 
3160  /** \brief Output a vector to a stream
3161 
3162  No trailing space is output after the last element, and an
3163  endline is output only if \c endline is set to \c true. If the
3164  parameter \c n is zero, this function silently does nothing.
3165 
3166  This works with any class <tt>vec_t</tt> which has an operator[]
3167  which returns either the value of or a reference to the ith
3168  element and the element type has its own output operator which
3169  has been defined.
3170  */
3171  template<class vec_t>
3172  void vector_out(std::ostream &os, const vec_t &v, bool endline=false) {
3173 
3174  size_t n=v.size();
3175 
3176  // This next line is important since n-1 is not well-defined if n=0
3177  if (n==0) return;
3178 
3179  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
3180  os << v[n-1];
3181  if (endline) os << std::endl;
3182  return;
3183  }
3184 
3185  /** \brief Fill a vector with a specified grid
3186  */
3187  template<class vec_t, class data_t>
3188  void vector_grid(uniform_grid<data_t> g, vec_t &v) {
3189  g.template vector<vec_t>(v);
3190  return;
3191  }
3192 
3193  /// Set a matrix to unity on the diagonal and zero otherwise
3194  template<class mat_t>
3195  void matrix_set_identity(size_t M, size_t N, mat_t &m) {
3196  for(size_t i=0;i<M;i++) {
3197  for(size_t j=0;j<N;j++) {
3198  if (i==j) m(i,j)=1.0;
3199  else m(i,j)=0.0;
3200  }
3201  }
3202  return;
3203  }
3204 
3205  /// Set a matrix to unity on the diagonal and zero otherwise
3206  template<class mat_t>
3207  void matrix_set_identity(mat_t &m) {
3208  matrix_set_identity(m.size1(),m.size2(),m);
3209  return;
3210  }
3211  //@}
3212 
3213  /// \name Vector range classes and functions in src/base/vector.h
3214  //@{
3215  /** \brief Vector range function for pointers
3216 
3217  \note In this case, the return type is the same as the
3218  type of the first parameter.
3219  */
3220  template<class dat_t> dat_t *vector_range
3221  (dat_t *v, size_t start, size_t last) {
3222  return v+start;
3223  }
3224 
3225  /** \brief Vector range function for const pointers
3226 
3227  \note In this case, the return type is the same as the
3228  type of the first parameter.
3229  */
3230  template<class dat_t> const dat_t *const_vector_range
3231  (const dat_t *v, size_t start, size_t last) {
3232  return v+start;
3233  }
3234 
3235  /** \brief Vector range function template for ublas vectors
3236 
3237  The element with index \c start in the original vector
3238  will become the first argument in the new vector, and
3239  the new vector will have size <tt>last-start</tt> .
3240 
3241  \note In this case, the return type is not the same as the
3242  type of the first parameter.
3243  */
3244  template<class dat_t> boost::numeric::ublas::vector_range
3247  size_t start, size_t last) {
3250  (v,boost::numeric::ublas::range(start,last));
3251  }
3252 
3253  /** \brief Const vector range function template for ublas vectors
3254 
3255  The element with index \c start in the original vector
3256  will become the first argument in the new vector, and
3257  the new vector will have size <tt>last-start</tt> .
3258 
3259  \note In this case, the return type is not the same as the
3260  type of the first parameter.
3261  */
3262  template<class dat_t> const boost::numeric::ublas::vector_range
3265  size_t start, size_t last) {
3268  (v,boost::numeric::ublas::range(start,last));
3269  }
3270 
3271  /** \brief Const vector range function template for const ublas
3272  vectors
3273 
3274  The element with index \c start in the original vector
3275  will become the first argument in the new vector, and
3276  the new vector will have size <tt>last-start</tt> .
3277 
3278  \note In this case, the return type is not the same as the
3279  type of the first parameter.
3280  */
3281  template<class dat_t> const boost::numeric::ublas::vector_range
3284  size_t start, size_t last) {
3287  (v,boost::numeric::ublas::range(start,last));
3288  }
3289 
3290  /** \brief Vector range function template for ublas vector
3291  ranges of ublas vectors
3292 
3293  The element with index \c start in the original vector
3294  will become the first argument in the new vector, and
3295  the new vector will have size <tt>last-start</tt> .
3296 
3297  \note In this case, the return type is not the same as the
3298  type of the first parameter.
3299  */
3300  template<class dat_t>
3304  vector_range
3307  size_t start, size_t last) {
3311  (v,boost::numeric::ublas::range(start,last));
3312  }
3313 
3314  /** \brief Const vector range function template for ublas vector
3315  ranges of ublas vectors
3316 
3317  The element with index \c start in the original vector
3318  will become the first argument in the new vector, and
3319  the new vector will have size <tt>last-start</tt> .
3320 
3321  \note In this case, the return type is not the same as the
3322  type of the first parameter.
3323  */
3324  template<class dat_t>
3331  size_t start, size_t last) {
3335  (v,boost::numeric::ublas::range(start,last));
3336  }
3337 
3338  /** \brief Const vector range function template for const
3339  ublas vector ranges of ublas vectors
3340 
3341  The element with index \c start in the original vector
3342  will become the first argument in the new vector, and
3343  the new vector will have size <tt>last-start</tt> .
3344 
3345  \note In this case, the return type is not the same as the
3346  type of the first parameter.
3347  */
3348  template<class dat_t>
3355  size_t start, size_t last) {
3359  (v,boost::numeric::ublas::range(start,last));
3360  }
3361 
3362  /** \brief Const vector range function template for const
3363  ublas vector ranges of const ublas vectors
3364 
3365  The element with index \c start in the original vector
3366  will become the first argument in the new vector, and
3367  the new vector will have size <tt>last-start</tt> .
3368 
3369  \note In this case, the return type is not the same as the
3370  type of the first parameter.
3371  */
3372  template<class dat_t>
3379  size_t start, size_t last) {
3383  (v,boost::numeric::ublas::range(start,last));
3384  }
3385 
3386  // Forward definition for friendship
3387  template<class vec_t> class const_vector_range_gen;
3388 
3389  /** \brief Experimental vector range object
3390  */
3391  template<class vec_t> class vector_range_gen {
3392 
3393  protected:
3394 
3395  friend class const_vector_range_gen<vec_t>;
3396 
3397  /// A reference to the original vector
3398  vec_t &v_;
3399 
3400  /// The index offset
3401  size_t start_;
3402 
3403  /// The end() iterator
3404  size_t last_;
3405 
3406  public:
3407 
3408  /// Create an object starting with index \c start in vector \c v
3409  vector_range_gen(vec_t &v, size_t start, size_t last) : v_(v),
3410  start_(start), last_(last) {
3411 #if !O2SCL_NO_RANGE_CHECK
3412  if (last<start) {
3413  O2SCL_ERR2("End before beginning in vector_range_gen::",
3414  "vector_range_gen(vec_t,size_t,size_t)",
3416  }
3417 #endif
3418  }
3419 
3420  /// Create an object from a previously constructed range object
3421  vector_range_gen(const vector_range_gen &v2, size_t start,
3422  size_t last) : v_(v2.v_),
3423  start_(start+v2.start_), last_(last+v2.start_) {
3424 #if !O2SCL_NO_RANGE_CHECK
3425  if (last<start) {
3426  O2SCL_ERR2("End before beginning in vector_range_gen::",
3427  "vector_range_gen(vector_range_gen,size_t,size_t)",
3429  }
3430  if (last>v2.last_) {
3431  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
3432  "vector_range_gen(vector_range_gen,size_t,size_t)",
3434  }
3435 #endif
3436  }
3437 
3438  /// Return the vector size
3439  size_t size() const {
3440  return last_-start_;
3441  }
3442 
3443  /// Return a reference ith element
3444  double &operator[](size_t i) {
3445 #if !O2SCL_NO_RANGE_CHECK
3446  if (i+start_>=last_) {
3447  O2SCL_ERR("Index out of range in vector_range_gen::operator[].",
3449  }
3450 #endif
3451  return v_[i+start_];
3452  }
3453 
3454  /// Return a const reference ith element
3455  const double &operator[](size_t i) const {
3456 #if !O2SCL_NO_RANGE_CHECK
3457  if (i+start_>=last_) {
3458  O2SCL_ERR2("Index out of range in ",
3459  "vector_range_gen::operator[] const.",o2scl::exc_einval);
3460  }
3461 #endif
3462  return v_[i+start_];
3463  }
3464  };
3465 
3466  /** \brief Experimental const vector range object
3467  */
3468  template<class vec_t> class const_vector_range_gen {
3469 
3470  protected:
3471 
3472  /// A reference to the original vector
3473  const vec_t &v_;
3474 
3475  /// The index offset
3476  size_t start_;
3477 
3478  /// The end() iterator
3479  size_t last_;
3480 
3481  public:
3482 
3483  /// Create an object starting with index \c start in vector \c v
3484  const_vector_range_gen(const vec_t &v, size_t start, size_t last) : v_(v),
3485  start_(start), last_(last) {
3486 #if !O2SCL_NO_RANGE_CHECK
3487  if (last<start) {
3488  O2SCL_ERR2("End before beginning in vector_range_gen::",
3489  "vector_range_gen(vec_t,size_t,size_t)",
3491  }
3492 #endif
3493  }
3494 
3495  /// Create an object from a previously constructed range object
3497  size_t last) : v_(v2.v_),
3498  start_(start+v2.start_), last_(last+v2.start_) {
3499 #if !O2SCL_NO_RANGE_CHECK
3500  if (last<start) {
3501  O2SCL_ERR2("End before beginning in vector_range_gen::",
3502  "vector_range_gen(vector_range_gen,size_t,size_t)",
3504  }
3505  if (last>v2.last_) {
3506  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
3507  "vector_range_gen(vector_range_gen,size_t,size_t)",
3509  }
3510 #endif
3511  }
3512 
3513  /// Create an object from a previously constructed range object
3515  size_t last) : v_(v2.v_),
3516  start_(start+v2.start_), last_(last+v2.start_) {
3517 #if !O2SCL_NO_RANGE_CHECK
3518  if (last<start) {
3519  O2SCL_ERR2("End before beginning in vector_range_gen::",
3520  "vector_range_gen(vector_range_gen,size_t,size_t)",
3522  }
3523  if (last>v2.last_) {
3524  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
3525  "vector_range_gen(vector_range_gen,size_t,size_t)",
3527  }
3528 #endif
3529  }
3530 
3531  /// Return the vector size
3532  size_t size() const {
3533  return last_-start_;
3534  }
3535 
3536  /// Return a const reference ith element
3537  const double &operator[](size_t i) const {
3538 #if !O2SCL_NO_RANGE_CHECK
3539  if (i+start_>=last_) {
3540  O2SCL_ERR2("Index out of range in ",
3541  "vector_range_gen::operator[] const.",o2scl::exc_einval);
3542  }
3543 #endif
3544  return v_[i+start_];
3545  }
3546  };
3547 
3548  /** \brief Create a \ref o2scl::vector_range_gen object
3549  from a <tt>std::vector</tt>
3550  */
3551  template<class data_t> vector_range_gen<std::vector<data_t> >
3552  vector_range(std::vector<data_t> &v, size_t start, size_t last) {
3553  return vector_range_gen<std::vector<data_t> >(v,start,last);
3554  }
3555 
3556  /** \brief Create a \ref o2scl::vector_range_gen object
3557  from a <tt>std::vector</tt>
3558  */
3559  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
3560  const_vector_range(const std::vector<data_t> &v, size_t start,
3561  size_t last) {
3562  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
3563  }
3564 
3565  /** \brief Create a \ref o2scl::vector_range_gen object
3566  from a <tt>std::vector</tt>
3567  */
3568  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
3569  const_vector_range(std::vector<data_t> &v, size_t start,
3570  size_t last) {
3571  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
3572  }
3573 
3574  /** \brief Recursively create a \ref o2scl::vector_range_gen object
3575  from a vector range
3576  */
3577  template<class vec_t> vector_range_gen<vec_t>
3578  vector_range(vector_range_gen<vec_t> &v, size_t start, size_t last) {
3579  return vector_range_gen<vec_t>(v,start,last);
3580  }
3581 
3582  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3583  object from a vector range
3584  */
3585  template<class vec_t> const const_vector_range_gen<vec_t>
3587  size_t start, size_t last) {
3588  return const_vector_range_gen<vec_t>(v,start,last);
3589  }
3590 
3591  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3592  object from a const vector range
3593  */
3594  template<class vec_t> const const_vector_range_gen<vec_t>
3596  size_t start, size_t last) {
3597  return const_vector_range_gen<vec_t>(v,start,last);
3598  }
3599 
3600  /** \brief Recursively create a const \ref o2scl::vector_range_gen
3601  object from a const vector range
3602  */
3603  template<class vec_t> const const_vector_range_gen<vec_t>
3605  size_t start, size_t last) {
3606  return const_vector_range_gen<vec_t>(v,start,last);
3607  }
3608 
3609  /** \brief Vector range function template for <tt>std::vector</tt>
3610 
3611  The element with index \c start in the original vector
3612  will become the first argument in the new vector, and
3613  the new vector will have size <tt>last-start</tt> .
3614 
3615  \note In this case, the return type is the same as the
3616  type of the first parameter.
3617  \note Unlike the ublas and pointer cases, this forces
3618  a copy.
3619  */
3620  template<class dat_t> std::vector<dat_t>
3621  vector_range_copy(const std::vector<dat_t> &v, size_t start, size_t last) {
3622  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3623  }
3624 
3625  /** \brief Const vector range function template for <tt>std::vector</tt>
3626 
3627  The element with index \c start in the original vector
3628  will become the first argument in the new vector, and
3629  the new vector will have size <tt>last-start</tt> .
3630 
3631  \note In this case, the return type is the same as the
3632  type of the first parameter.
3633  \note Unlike the ublas and pointer cases, this forces
3634  a copy.
3635  */
3636  template<class dat_t> const std::vector<dat_t>
3637  vector_range_copy(const std::vector<dat_t> &v, size_t start,
3638  size_t last) {
3639  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3640  }
3641 
3642  //@}
3643 
3644 }
3645 
3646 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
3647 
3648 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
3649 #include <armadillo>
3650 namespace o2scl {
3651 
3652  /// \name Armadillo specializations in src/base/vector.h
3653  //@{
3654  /// Armadillo version of \ref matrix_max()
3655  double matrix_max(const arma::mat &data);
3656 
3657  /// Armadillo version of \ref matrix_min()
3658  double matrix_min(const arma::mat &data);
3659 
3660  /// Armadillo version of \ref matrix_row()
3661  template<> arma::subview_row<double>
3662  matrix_row<arma::mat,arma::subview_row<double> >
3663  (arma::mat &M, size_t row);
3664 
3665  /// Armadillo version of \ref matrix_column()
3666  template<> arma::subview_col<double>
3667  matrix_column<arma::mat,arma::subview_col<double> >
3668  (arma::mat &M, size_t column);
3669  //@}
3670 
3671 }
3672 
3673 #endif
3674 
3675 #if defined (O2SCL_EIGEN) || defined (DOXYGEN)
3676 #include <eigen3/Eigen/Dense>
3677 
3678 namespace o2scl {
3679 
3680  /// \name Eigen specializations in src/base/vector.h
3681  //@{
3682  /// Eigen version of \ref matrix_max()
3683  double matrix_max(const Eigen::MatrixXd &data);
3684 
3685  /// Eigen version of \ref matrix_min()
3686  double matrix_min(const Eigen::MatrixXd &data);
3687 
3688  /// Eigen version of \ref matrix_row()
3689  template<> Eigen::MatrixXd::RowXpr
3690  matrix_row<Eigen::MatrixXd,Eigen::MatrixXd::RowXpr>
3691  (Eigen::MatrixXd &M, size_t row);
3692 
3693  /// Eigen version of \ref matrix_column()
3694  template<> Eigen::MatrixXd::ColXpr
3695  matrix_column<Eigen::MatrixXd,Eigen::MatrixXd::ColXpr>
3696  (Eigen::MatrixXd &M, size_t column);
3697  //@}
3698 
3699 }
3700 
3701 #endif
3702 
3703 #else
3704 
3705 #include <o2scl/vector_special.h>
3706 
3707 // End of "#if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)"
3708 #endif
3709 
3710 #ifdef DOXYGEN
3711 /** \brief Placeholder documentation of some related Boost objects
3712  */
3713 namespace boost {
3714  /** \brief Documentation of Boost::numeric objects
3715  */
3716  namespace numeric {
3717  /** \brief Documentation of uBlas objects
3718  */
3719  namespace ublas {
3720  /** \brief The default vector type from uBlas
3721 
3722  The uBlas types aren't documented here, but the full documentation
3723  is available at
3724  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3725 
3726  Internally in \o2, this is often typedef'd using
3727  \code
3728  typedef boost::numeric::ublas::vector<double> ubvector;
3729  typedef boost::numeric::ublas::vector<size_t> ubvector_size_t;
3730  typedef boost::numeric::ublas::vector<int> ubvector_int;
3731  \endcode
3732 
3733  This is documented in \ref vector.h .
3734  */
3735  template<class T, class A> class vector {
3736  };
3737  /** \brief The default matrix type from uBlas
3738 
3739  The uBlas types aren't documented here, but the full documentation
3740  is available at
3741  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3742 
3743  Internally in \o2, this is often typedef'd using
3744  \code
3745  typedef boost::numeric::ublas::matrix<double> ubmatrix;
3746  typedef boost::numeric::ublas::matrix<size_t> ubmatrix_size_t;
3747  typedef boost::numeric::ublas::matrix<int> ubmatrix_int;
3748  \endcode
3749 
3750  This is documented in \ref vector.h .
3751  */
3752  template<class T, class F, class A> class matrix {
3753  };
3754  }
3755  }
3756 }
3757 // End of "#ifdef DOXYGEN"
3758 #endif
3759 
3760 // End of "#ifndef O2SCL_VECTOR_H"
3761 #endif
o2scl::matrix_view_omit_row::matrix_view_omit_row
matrix_view_omit_row(mat_t &m, size_t row_omit)
Create.
Definition: vector.h:2768
o2scl::matrix_view_omit_row::operator()
const double & operator()(size_t i, size_t j) const
Return a const reference.
Definition: vector.h:2781
boost::numeric::ublas::matrix
The default matrix type from uBlas.
Definition: vector.h:3752
o2scl::vector_norm_double
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.
Definition: vector.h:2296
o2scl::vector_range_gen::last_
size_t last_
The end() iterator.
Definition: vector.h:3404
o2scl::matrix_view_vec_vec::swap
friend void swap(matrix_view_vec_vec &t1, matrix_view_vec_vec &t2)
Swap method.
Definition: vector.h:2953
o2scl::matrix_row_gen_ctor::matrix_row_gen_ctor
matrix_row_gen_ctor(size_t n_cols=0)
Create a row object from row row of matrix m.
Definition: vector.h:2658
o2scl::matrix_row_gen::operator[]
double & operator[](size_t i)
Return a reference to the ith column of the selected row.
Definition: vector.h:2624
o2scl::vector_bsearch
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.
Definition: vector.h:1977
o2scl::matrix_row_gen_ctor
Matrix row object with a constructor and resize method.
Definition: vector.h:2638
o2scl::matrix_view_transpose::m_
mat_t & m_
A reference to the original matrix.
Definition: vector.h:2719
o2scl::matrix_swap_rows
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.
Definition: vector.h:623
o2scl::vector_range_gen::v_
vec_t & v_
A reference to the original vector.
Definition: vector.h:3398
o2scl::vector_index_vector
Trivial index vector.
Definition: vector.h:2525
o2scl::matrix_max_value
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.
Definition: vector.h:1425
boost::numeric::ublas::vector
The default vector type from uBlas.
Definition: vector.h:3735
o2scl::const_matrix_column_gen::column_
size_t column_
The selected column.
Definition: vector.h:3116
o2scl::const_vector_range_gen::size
size_t size() const
Return the vector size.
Definition: vector.h:3532
o2scl::vector_index_vector_size::size
size_t size() const
Get the size of the vector.
Definition: vector.h:2564
o2scl::vector_range_gen::size
size_t size() const
Return the vector size.
Definition: vector.h:3439
o2scl::matrix_view_omit_row::size1
size_t size1() const
Return the number of rows.
Definition: vector.h:2790
o2scl::matrix_view_omit_row
Construct a view of a matrix omtting a specified row.
Definition: vector.h:2756
o2scl::const_vector_range_gen::start_
size_t start_
The index offset.
Definition: vector.h:3476
o2scl::matrix_copy
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
Definition: vector.h:176
o2scl::matrix_minmax_index
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.
Definition: vector.h:1712
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::const_matrix_view::operator()
const double & operator()(size_t row, size_t col) const
Return a reference to the element at row row and column col.
o2scl::matrix_view_omit_row::operator()
double & operator()(size_t i, size_t j)
Return a reference.
Definition: vector.h:2773
o2scl::matrix_transpose
void matrix_transpose(mat_t &src, mat2_t &dest)
Simple transpose.
Definition: vector.h:222
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::matrix_is_lower
bool matrix_is_lower(mat_t &src)
Simple test that a matrix is lower triangular.
Definition: vector.h:303
o2scl::vector_minmax
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.
Definition: vector.h:1322
o2scl::matrix_min_index
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.
Definition: vector.h:1622
o2scl::vector_range_gen::operator[]
double & operator[](size_t i)
Return a reference ith element.
Definition: vector.h:3444
o2scl::matrix_swap_rows_double
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.
Definition: vector.h:640
o2scl::matrix_max
double matrix_max(const arma::mat &data)
Armadillo version of matrix_max()
o2scl::const_matrix_row_gen
Generic object which represents a row of a const matrix.
Definition: vector.h:2864
o2scl::matrix_row_gen::m_
mat_t & m_
A reference to the original matrix.
Definition: vector.h:2612
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::matrix_swap_double
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.
Definition: vector.h:561
o2scl::vector_range
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.
Definition: vector.h:3578
o2scl::vector_index_vector_size::vector_index_vector_size
vector_index_vector_size(size_t n_)
Create an index vector with size n_.
Definition: vector.h:2549
o2scl::gsl_vector_wrap
A simple convenience wrapper for GSL vector objects.
Definition: vector.h:68
o2scl::vector_min_value
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1196
o2scl::vector_lookup
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.
Definition: vector.h:1780
o2scl::matrix_view
A simple matrix view object.
Definition: vector.h:2909
o2scl::matrix_column
mat_column_t matrix_column(mat_t &M, size_t column)
Construct a column of a matrix.
Definition: vector.h:3051
o2scl::vector_max_value
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
Definition: vector.h:1121
boost
Placeholder documentation of some related Boost objects.
Definition: vector.h:3713
o2scl::const_matrix_column_gen::operator[]
const double & operator[](size_t i) const
Return a const reference to the ith row of the selected column.
Definition: vector.h:3126
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::const_matrix_row_gen::const_matrix_row_gen
const_matrix_row_gen(const mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2877
o2scl::const_matrix_column_gen::m_
const mat_t & m_
A reference to the original matrix.
Definition: vector.h:3113
o2scl::matrix_view_omit_column::operator()
double & operator()(size_t i, size_t j)
Return a reference.
Definition: vector.h:2824
o2scl::vector_swap_double
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.
Definition: vector.h:490
o2scl::vector_min_quad_loc
data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector minimum by quadratic fit.
Definition: vector.h:1412
o2scl::matrix_view_transpose::operator()
double & operator()(size_t i, size_t j)
Return a reference to the ith column of the selected row.
Definition: vector.h:2728
o2scl::matrix_column_gen::matrix_column_gen
matrix_column_gen(mat_t &m, size_t column)
Create a column object from column column of matrix m.
Definition: vector.h:3082
o2scl::matrix_min
double matrix_min(const arma::mat &data)
Armadillo version of matrix_min()
o2scl::matrix_view::size1
size_t size1() const
Return the number of rows.
o2scl::const_matrix_column_gen::const_matrix_column_gen
const_matrix_column_gen(const mat_t &m, size_t column)
Create a column object from column column of matrix m.
Definition: vector.h:3121
o2scl::matrix_view_omit_row::size2
size_t size2() const
Return the number of columns.
Definition: vector.h:2796
o2scl::const_vector_range_gen::const_vector_range_gen
const_vector_range_gen(const vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
Definition: vector.h:3484
o2scl::vector_set_all
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.
Definition: vector.h:2338
o2scl::vector_minmax_value
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.
Definition: vector.h:1272
o2scl::const_matrix_view::size1
size_t size1() const
Return the number of rows.
o2scl::matrix_row_gen_ctor::matrix_row_gen_ctor
matrix_row_gen_ctor(mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2654
o2scl::vector_range_gen::vector_range_gen
vector_range_gen(vec_t &v, size_t start, size_t last)
Create an object starting with index start in vector v.
Definition: vector.h:3409
o2scl::matrix_view_vec_vec
View a o2scl::table object as a matrix.
Definition: vector.h:2942
o2scl::matrix_view_vec_vec::matrix_view_vec_vec
matrix_view_vec_vec(vec2_t &vv)
Create a matrix view object from the specified table and list of rows.
Definition: vector.h:2967
o2scl::vector_range_gen::vector_range_gen
vector_range_gen(const vector_range_gen &v2, size_t start, size_t last)
Create an object from a previously constructed range object.
Definition: vector.h:3421
o2scl::matrix_row_gen_ctor::row_
size_t row_
The selected row.
Definition: vector.h:2646
o2scl::matrix_row_gen_ctor::mp
mat_t * mp
A pointer to the matrix.
Definition: vector.h:2643
o2scl::matrix_column_gen::column_
size_t column_
The selected column.
Definition: vector.h:3077
o2scl::matrix_make_upper
void matrix_make_upper(mat_t &src)
Make a matrix upper triangular by setting the lower triangular entries to zero.
Definition: vector.h:346
o2scl::vector_sort_index
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)
Definition: vector.h:800
o2scl::matrix_view_omit_column
Construct a view of a matrix omitting one specified column.
Definition: vector.h:2807
o2scl::matrix_minmax
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.
Definition: vector.h:1679
o2scl::matrix_row
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
Definition: vector.h:2593
o2scl::const_matrix_view::size2
size_t size2() const
Return the number of columns.
o2scl::matrix_sum
data_t matrix_sum(size_t m, size_t n, const mat_t &data)
Compute the sum of matrix elements.
Definition: vector.h:1745
o2scl::const_vector_range_gen::operator[]
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:3537
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::vector_grid
void vector_grid(uniform_grid< data_t > g, vec_t &v)
Fill a vector with a specified grid.
Definition: vector.h:3188
o2scl::matrix_view_omit_column::matrix_view_omit_column
matrix_view_omit_column(mat_t &m, size_t column_omit)
Create.
Definition: vector.h:2819
o2scl::const_matrix_row_gen::operator[]
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2881
o2scl::vector_range_gen
Experimental vector range object.
Definition: vector.h:3391
o2scl::vector_norm
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...
Definition: vector.h:2252
o2scl::vector_bsearch_inc
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.
Definition: vector.h:1906
o2scl::matrix_min_value_double
double matrix_min_value_double(const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1597
o2scl::sort_downheap
void sort_downheap(vec_t &data, size_t n, size_t k)
Provide a downheap() function for vector_sort()
Definition: vector.h:650
o2scl::matrix_row_gen
Generic object which represents a row of a matrix.
Definition: vector.h:2607
o2scl::matrix_is_upper
bool matrix_is_upper(mat_t &src)
Simple test that a matrix is upper triangular.
Definition: vector.h:317
o2scl::matrix_view_vec_vec::vvp
vec2_t * vvp
Pointer to the table.
Definition: vector.h:2947
o2scl::vector_max_quad_loc
data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y)
Location of vector maximum by quadratic fit.
Definition: vector.h:1377
o2scl::vector_sum_double
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.
Definition: vector.h:2223
o2scl::matrix_column_gen::operator[]
const double & operator[](size_t i) const
Return a const reference to the ith row of the selected column.
Definition: vector.h:3091
o2scl::const_matrix_row_gen::row_
size_t row_
The selected row.
Definition: vector.h:2872
o2scl::vector_index_vector_size::resize
void resize(size_t n_)
Resize the index vector.
Definition: vector.h:2570
o2scl::matrix_view_vec_vec::operator()
double & operator()(size_t row, size_t col)
Return a reference to the element at row row and column col.
Definition: vector.h:3011
o2scl::vector_sort
void vector_sort(size_t n, vec_t &data)
Sort a vector (in increasing order)
Definition: vector.h:705
o2scl::const_matrix_column_gen
Generic object which represents a column of a const matrix.
Definition: vector.h:3108
o2scl::matrix_view_vec_vec::size2
size_t size2() const
Return the number of columns.
Definition: vector.h:2980
o2scl::vector_max_quad
data_t vector_max_quad(size_t n, const vec_t &data)
Maximum of vector by quadratic fit.
Definition: vector.h:1352
o2scl::matrix_view_vec_vec::operator()
const double & operator()(size_t row, size_t col) const
Return a reference to the element at row row and column col.
Definition: vector.h:2989
o2scl::vector_index_vector_size::operator[]
data_t operator[](size_t &i) const
Obtain the element with index i.
Definition: vector.h:2555
o2scl::vector_index_vector_size::n
size_t n
The vector size.
Definition: vector.h:2543
o2scl::matrix_view_transpose::matrix_view_transpose
matrix_view_transpose(mat_t &m)
Create a row object from row row of matrix m.
Definition: vector.h:2724
o2scl::const_vector_range
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
Definition: vector.h:3231
o2scl::matrix_column_gen::m_
mat_t & m_
A reference to the original matrix.
Definition: vector.h:3074
o2scl::matrix_make_lower
void matrix_make_lower(mat_t &src)
Make a matrix lower triangular by setting the upper triangular entries to zero.
Definition: vector.h:332
o2scl::vector_reverse_double
void vector_reverse_double(size_t n, vec_t &data)
Reverse the first n elements in a vector of double precision numbers.
Definition: vector.h:2491
o2scl::matrix_min_value
data_t matrix_min_value(size_t m, size_t n, const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1551
o2scl::uniform_grid
A class representing a uniform linear or logarithmic grid.
Definition: uniform_grid.h:38
o2scl::matrix_row_gen_ctor::operator[]
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2701
o2scl::vector_is_monotonic
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.
Definition: vector.h:2019
o2scl::vector_rotate
void vector_rotate(size_t n, vec_t &data, size_t k)
"Rotate" a vector so that the kth element is now the beginning
Definition: vector.h:2437
o2scl::sort_index_downheap
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()
Definition: vector.h:734
o2scl::matrix_view_omit_row::m_
mat_t & m_
A reference to the original matrix.
Definition: vector.h:2761
o2scl::vector_is_finite
bool vector_is_finite(size_t n, vec_t &data)
Test if the first n elements of a vector are finite.
Definition: vector.h:2152
o2scl::vector_copy_jackknife
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.
Definition: vector.h:2379
o2scl::vector_range_copy
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
Definition: vector.h:3621
o2scl::vector_is_strictly_monotonic
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 ...
Definition: vector.h:2101
o2scl::matrix_swap_cols_double
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.
Definition: vector.h:612
o2scl::matrix_view::operator()
const double & operator()(size_t row, size_t col) const
Return a reference to the element at row row and column col.
o2scl::const_matrix_row_gen::m_
const mat_t & m_
A reference to the original matrix.
Definition: vector.h:2869
o2scl::const_vector_range_gen::const_vector_range_gen
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.
Definition: vector.h:3514
o2scl::vector_min
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.
Definition: vector.h:1251
o2scl::itos
std::string itos(int x)
Convert an integer to a string.
o2scl::vector_max
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.
Definition: vector.h:1176
o2scl::matrix_view_transpose
Construct a view of the transpose of a matrix.
Definition: vector.h:2714
o2scl::matrix_lookup
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.
Definition: vector.h:1827
o2scl::vector_smallest
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.
Definition: vector.h:920
o2scl::matrix_row_gen_ctor::resize
void resize(size_t n_cols=0)
Resize.
Definition: vector.h:2670
o2scl::vector_minmax_index
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.
Definition: vector.h:1295
o2scl::matrix_view_omit_column::m_
mat_t & m_
A reference to the original matrix.
Definition: vector.h:2812
o2scl::gsl_matrix_wrap
A simple convenience wrapper for GSL matrix objects.
Definition: vector.h:89
o2scl::matrix_set_all
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.
Definition: vector.h:2356
o2scl::vector_index_vector_size
Index vector with a size method.
Definition: vector.h:2538
o2scl::vector_sort_double
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Definition: vector.h:895
o2scl::vector_smallest_index
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.
Definition: vector.h:993
o2scl::const_vector_range_gen
Experimental const vector range object.
Definition: vector.h:3387
o2scl::vector_range
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:3221
o2scl::matrix_row_gen::row_
size_t row_
The selected row.
Definition: vector.h:2615
o2scl::matrix_view_transpose::size2
size_t size2() const
Return the number of columns.
Definition: vector.h:2745
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::matrix_column_gen
Generic object which represents a column of a matrix.
Definition: vector.h:3069
o2scl::const_vector_range_gen::const_vector_range_gen
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.
Definition: vector.h:3496
o2scl::matrix_set_identity
void matrix_set_identity(size_t M, size_t N, mat_t &m)
Set a matrix to unity on the diagonal and zero otherwise.
Definition: vector.h:3195
o2scl::vector_sum
data_t vector_sum(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector.
Definition: vector.h:2180
o2scl::vector_min_quad
data_t vector_min_quad(size_t n, const vec_t &data)
Minimum of vector by quadratic fit.
Definition: vector.h:1387
o2scl::vector_largest
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.
Definition: vector.h:1058
o2scl::vector_reverse
void vector_reverse(size_t n, vec_t &data)
Reverse the first n elements of a vector.
Definition: vector.h:2456
o2scl::matrix_view_transpose::size1
size_t size1() const
Return the number of rows.
Definition: vector.h:2739
o2scl::const_matrix_view
A simple matrix view object.
Definition: vector.h:2888
o2scl::matrix_view_transpose::operator()
const double & operator()(size_t i, size_t j) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2733
o2scl::vector_range_gen::start_
size_t start_
The index offset.
Definition: vector.h:3401
o2scl::matrix_max_value_double
double matrix_max_value_double(const mat_t &data)
Compute the maximum of a matrix.
Definition: vector.h:1470
o2scl::matrix_view_vec_vec::size1
size_t size1() const
Return the number of rows.
Definition: vector.h:2973
o2scl::vector_diffs
void vector_diffs(const vec_t &v_data, rvec_t &v_diffs)
Create a new vector containing the differences between adjacent entries.
Definition: vector.h:2193
o2scl::matrix_row_gen_ctor::size
size_t size() const
Return size.
Definition: vector.h:2684
o2scl::matrix_view_omit_column::size1
size_t size1() const
Return the number of rows.
Definition: vector.h:2841
o2scl::matrix_swap
void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2)
Swap of the first elements in two matrices.
Definition: vector.h:542
o2scl::const_vector_range_gen::v_
const vec_t & v_
A reference to the original vector.
Definition: vector.h:3473
o2scl::matrix_max_index
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.
Definition: vector.h:1494
o2scl::vector_bsearch_dec
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.
Definition: vector.h:1951
o2scl::matrix_view::size2
size_t size2() const
Return the number of columns.
o2scl::matrix_row_gen::operator[]
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2629
o2scl::matrix_row_gen_ctor::operator[]
double & operator[](size_t i)
Return a reference to the ith column of the selected row.
Definition: vector.h:2692
o2scl::matrix_swap_cols
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.
Definition: vector.h:595
o2scl::vector_out
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.
Definition: vector.h:3145
o2scl::matrix_view_omit_column::operator()
const double & operator()(size_t i, size_t j) const
Return a const reference.
Definition: vector.h:2832
o2scl::matrix_column_gen::operator[]
double & operator[](size_t i)
Return a reference to the ith row of the selected column.
Definition: vector.h:3086
o2scl::matrix_row_gen::matrix_row_gen
matrix_row_gen(mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2620
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::matrix_row_gen_ctor::mat
mat_t mat
A matrix to point to.
Definition: vector.h:2649
o2scl::vector_range_gen::operator[]
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:3455
o2scl::vector_max_index
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.
Definition: vector.h:1157
o2scl::matrix_view_omit_column::size2
size_t size2() const
Return the number of columns.
Definition: vector.h:2847
o2scl::vector_swap
void vector_swap(size_t N, vec_t &v1, vec2_t &v2)
Swap the first N elements of two vectors.
Definition: vector.h:420
o2scl::vector_min_index
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.
Definition: vector.h:1232
o2scl::const_vector_range_gen::last_
size_t last_
The end() iterator.
Definition: vector.h:3479

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).