vector.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, 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  \future Create a matrix transpose copy function?
49  \future Create matrix swap row and column functions
50 */
51 #include <iostream>
52 #include <cmath>
53 #include <string>
54 #include <fstream>
55 #include <sstream>
56 
57 #include <gsl/gsl_vector.h>
58 #include <gsl/gsl_sys.h>
59 #include <gsl/gsl_matrix.h>
60 
61 #include <o2scl/misc.h>
62 #include <o2scl/uniform_grid.h>
63 
64 namespace o2scl {
65 
66  /** \brief A simple convenience wrapper for GSL vector objects
67 
68  \warning This uses typecasts on externally allocated GSL
69  pointers and is not safe or fully const-correct.
70  */
72  const double *d;
73  size_t sz;
74  public:
75  gsl_vector_wrap(gsl_vector *m) {
76  d=(const double *)m->data;
77  sz=m->size;
78  }
79  double operator[](size_t i) const {
80  return d[i];
81  }
82  size_t size() {
83  return sz;
84  }
85  };
86 
87  /** \brief A simple convenience wrapper for GSL matrix objects
88 
89  \warning This uses typecasts on externally allocated GSL
90  pointers and is not safe or fully const-correct.
91  */
93  const double *d;
94  size_t sz1;
95  size_t sz2;
96  size_t tda;
97  public:
98  gsl_matrix_wrap(gsl_matrix *m) {
99  d=(const double *)m->data;
100  sz1=m->size1;
101  sz2=m->size2;
102  tda=m->tda;
103  }
104  double operator()(size_t i, size_t j) const {
105  return *(d+i*tda+j);
106  }
107  size_t size1() {
108  return sz1;
109  }
110  size_t size2() {
111  return sz2;
112  }
113  };
114 
115  /// \name Copying vectors and matrices
116  //@{
117  /** \brief Simple vector copy
118 
119  Copy \c src to \c dest, resizing \c dest if it is too small
120  to hold <tt>src.size()</tt> elements.
121 
122  This function will work for any classes \c vec_t and
123  \c vec2_t which have suitably defined <tt>operator[]</tt>,
124  <tt>size()</tt>, and <tt>resize()</tt> methods.
125  */
126  template<class vec_t, class vec2_t>
127  void vector_copy(const vec_t &src, vec2_t &dest) {
128  size_t N=src.size();
129  if (dest.size()<N) dest.resize(N);
130  size_t i, m=N%4;
131  for(i=0;i<m;i++) {
132  dest[i]=src[i];
133  }
134  for(i=m;i+3<N;i+=4) {
135  dest[i]=src[i];
136  dest[i+1]=src[i+1];
137  dest[i+2]=src[i+2];
138  dest[i+3]=src[i+3];
139  }
140  return;
141  }
142 
143  /** \brief Simple vector copy of the first N elements
144 
145  Copy the first \c N elements of \c src to \c dest.
146  It is assumed that the memory allocation for \c dest
147  has already been performed.
148 
149  This function will work for any class <tt>vec2_t</tt> which has
150  an operator[] which returns a reference to the corresponding
151  element and class <tt>vec_t</tt> with an operator[] which
152  returns either a reference or the value of the corresponding
153  element.
154  */
155  template<class vec_t, class vec2_t>
156  void vector_copy(size_t N, const vec_t &src, vec2_t &dest) {
157  size_t i, m=N%4;
158  for(i=0;i<m;i++) {
159  dest[i]=src[i];
160  }
161  for(i=m;i+3<N;i+=4) {
162  dest[i]=src[i];
163  dest[i+1]=src[i+1];
164  dest[i+2]=src[i+2];
165  dest[i+3]=src[i+3];
166  }
167  return;
168  }
169 
170  /** \brief Simple matrix copy
171 
172  Copy \c src to \c dest, resizing \c dest if it is too small.
173 
174  This function will work for any classes \c mat_t and
175  \c mat2_t which have suitably defined <tt>operator()</tt>,
176  <tt>size()</tt>, and <tt>resize()</tt> methods.
177  */
178  template<class mat_t, class mat2_t>
179  void matrix_copy(mat_t &src, mat2_t &dest) {
180  size_t m=src.size1();
181  size_t n=src.size2();
182  if (dest.size1()<m || dest.size2()<n) dest.resize(m,n);
183  for(size_t i=0;i<m;i++) {
184  for(size_t j=0;j<n;j++) {
185  dest(i,j)=src(i,j);
186  }
187  }
188  }
189 
190  /** \brief Simple matrix copy of the first \f$ (M,N) \f$
191  matrix elements
192 
193  Copy the first <tt>(M,N)</tt> elements of \c src to \c dest. It
194  is assumed that the memory allocation for \c dest has already
195  been performed.
196 
197  This function will work for any class <tt>vec2_t</tt> which has
198  an operator[][] which returns a reference to the corresponding
199  element and class <tt>vec_t</tt> with an operator[][] which
200  returns either a reference or the value of the corresponding
201  element.
202  */
203  template<class mat_t, class mat2_t>
204  void matrix_copy(size_t M, size_t N, mat_t &src, mat2_t &dest) {
205  for(size_t i=0;i<M;i++) {
206  for(size_t j=0;j<N;j++) {
207  dest(i,j)=src(i,j);
208  }
209  }
210  }
211  //@}
212 
213  /// \name Tranpositions
214  //@{
215  /** \brief Simple transpose
216 
217  Copy the transpose of \c src to \c dest, resizing \c dest if it
218  is too small.
219 
220  This function will work for any classes \c mat_t and
221  \c mat2_t which have suitably defined <tt>operator()</tt>,
222  <tt>size()</tt>, and <tt>resize()</tt> methods.
223  */
224  template<class mat_t, class mat2_t>
225  void matrix_transpose(mat_t &src, mat2_t &dest) {
226  size_t m=src.size1();
227  size_t n=src.size2();
228  if (dest.size1()<n || dest.size2()<m) dest.resize(n,m);
229  for(size_t i=0;i<m;i++) {
230  for(size_t j=0;j<n;j++) {
231  dest(i,j)=src(j,i);
232  }
233  }
234  }
235 
236  /** \brief Simple transpose of the first \f$ (m,n) \f$
237  matrix elements
238 
239  Copy the transpose of the first \c m rows and the first \c cols
240  of the matrix \c src into the matrix \c dest
241 
242  This function will work for any classes \c mat_t and \c mat2_t
243  which has a suitably defined <tt>operator()</tt> method.
244  */
245  template<class mat_t, class mat2_t>
246  void matrix_transpose(size_t m, size_t n, mat_t &src, mat2_t &dest) {
247  for(size_t i=0;i<m;i++) {
248  for(size_t j=0;j<n;j++) {
249  dest(i,j)=src(j,i);
250  }
251  }
252  }
253 
254  /** \brief Simple transpose in-place
255 
256  Transpose the matrix \c src . If the matrix is not square,
257  only the upper-left square part of the matrix will be
258  transposed.
259 
260  This function will work for any classes \c mat_t and
261  \c mat2_t which have suitably defined <tt>operator()</tt>,
262  <tt>size()</tt>, and <tt>resize()</tt> methods.
263  */
264  template<class mat_t, class data_t>
265  void matrix_transpose(mat_t &src) {
266  size_t m=src.size1();
267  size_t n=src.size2();
268  // Choose the smaller of n and m
269  if (m<n) n=m;
270  data_t tmp;
271  for(size_t i=0;i<n;i++) {
272  for(size_t j=0;j<n;j++) {
273  tmp=src(i,j);
274  src(i,j)=src(j,i);
275  src(j,i)=tmp;
276  }
277  }
278  }
279 
280  /** \brief Simple in-place transpose of the first \f$ (m,n) \f$
281  matrix elements
282 
283  Copy the transpose of the first \c m rows and the first \c cols
284  of the matrix \c src into the matrix \c dest
285 
286  This function will work for any classes \c mat_t and \c mat2_t
287  which has a suitably defined <tt>operator()</tt> method.
288  */
289  template<class mat_t, class data_t>
290  void matrix_transpose(size_t m, size_t n, mat_t &src) {
291  data_t tmp;
292  for(size_t i=0;i<m;i++) {
293  for(size_t j=0;j<n;j++) {
294  tmp=src(i,j);
295  src(i,j)=src(j,i);
296  src(j,i)=tmp;
297  }
298  }
299  }
300  //@}
301 
302  /// \name Upper and lower triangular functions
303  //@{
304  /** \brief Simple test that a matrix is lower triangular
305  */
306  template<class mat_t> bool matrix_is_lower(mat_t &src) {
307  size_t m=src.size1();
308  size_t n=src.size2();
309  bool ret=true;
310  for(size_t i=0;i<m;i++) {
311  for(size_t j=i+1;j<n;j++) {
312  if (src(i,j)!=0.0) return false;
313  }
314  }
315  return ret;
316  }
317 
318  /** \brief Simple test that a matrix is upper triangular
319  */
320  template<class mat_t> bool matrix_is_upper(mat_t &src) {
321  size_t m=src.size1();
322  size_t n=src.size2();
323  bool ret=true;
324  for(size_t j=0;j<n;j++) {
325  for(size_t i=j+1;i<m;i++) {
326  if (src(i,j)!=0.0) return false;
327  }
328  }
329  return ret;
330  }
331 
332  /** \brief Make a matrix lower triangular by setting the upper
333  triangular entries to zero
334  */
335  template<class mat_t> void matrix_make_lower(mat_t &src) {
336  size_t m=src.size1();
337  size_t n=src.size2();
338  for(size_t i=0;i<m;i++) {
339  for(size_t j=i+1;j<n;j++) {
340  src(i,j)=0.0;
341  }
342  }
343  return;
344  }
345 
346  /** \brief Make a matrix upper triangular by setting the lower
347  triangular entries to zero
348  */
349  template<class mat_t> void matrix_make_upper(mat_t &src) {
350  size_t m=src.size1();
351  size_t n=src.size2();
352  for(size_t j=0;j<n;j++) {
353  for(size_t i=j+1;i<m;i++) {
354  src(i,j)=0.0;
355  }
356  }
357  return;
358  }
359 
360  /** \brief Simple test that a matrix is lower triangular
361  for the first \c m rows and \c n columns
362  */
363  template<class mat_t> bool matrix_is_lower(size_t m, size_t n,
364  mat_t &src) {
365  bool ret=true;
366  for(size_t i=0;i<m;i++) {
367  for(size_t j=i+1;j<n;j++) {
368  if (src(i,j)!=0.0) return false;
369  }
370  }
371  return ret;
372  }
373 
374  /** \brief Simple test that a matrix is upper triangular
375  for the first \c m rows and \c n columns
376  */
377  template<class mat_t> bool matrix_is_upper(size_t m, size_t n,
378  mat_t &src) {
379  bool ret=true;
380  for(size_t j=0;j<n;j++) {
381  for(size_t i=j+1;i<m;i++) {
382  if (src(i,j)!=0.0) return false;
383  }
384  }
385  return ret;
386  }
387 
388  /** \brief Make the first \c m rows and \c n columns of a matrix
389  lower triangular by setting the upper triangular entries to zero
390  */
391  template<class mat_t> void matrix_make_lower(size_t m, size_t n,
392  mat_t &src) {
393  for(size_t i=0;i<m;i++) {
394  for(size_t j=i+1;j<n;j++) {
395  src(i,j)=0.0;
396  }
397  }
398  return;
399  }
400 
401  /** \brief Make the first \c m rows and \c n columns of a matrix
402  upper triangular by setting the lower triangular entries to zero
403  */
404  template<class mat_t> void matrix_make_upper(size_t m, size_t n,
405  mat_t &src) {
406  for(size_t j=0;j<n;j++) {
407  for(size_t i=j+1;i<m;i++) {
408  src(i,j)=0.0;
409  }
410  }
411  return;
412  }
413  //@}
414 
415  /// \name Swapping parts of vectors and matrices
416  //@{
417  /** \brief Swap the first \c N elements of two vectors
418 
419  This function swaps the elements of \c v1 and \c v2, one element
420  at a time.
421  */
422  template<class vec_t, class vec2_t, class data_t>
423  void vector_swap(size_t N, vec_t &v1, vec2_t &v2) {
424  data_t temp;
425  size_t i, m=N%4;
426  for(i=0;i<m;i++) {
427  temp=v1[i];
428  v1[i]=v2[i];
429  v2[i]=temp;
430  }
431  for(i=m;i+3<N;i+=4) {
432  temp=v1[i];
433  v1[i]=v2[i];
434  v2[i]=temp;
435  temp=v1[i+1];
436  v1[i+1]=v2[i+1];
437  v2[i+1]=temp;
438  temp=v1[i+2];
439  v1[i+2]=v2[i+2];
440  v2[i+2]=temp;
441  temp=v1[i+3];
442  v1[i+3]=v2[i+3];
443  v2[i+3]=temp;
444  }
445  return;
446  }
447 
448  /** \brief Swap all elements in two vectors
449 
450  This function swaps the elements of \c v1 and \c v2, one element
451  at a time.
452 
453  \note It is almost always better to use <tt>std::swap</tt>
454  than this function, which is provided only in cases where
455  one knows one is going to be forced to use a vector type
456  without a properly defined <tt>std::swap</tt> method.
457  */
458  template<class vec_t, class vec2_t, class data_t>
459  void vector_swap(vec_t &v1, vec2_t &v2) {
460  size_t N=v1.size();
461  if (v2.size()<N) N=v2.size();
462  data_t temp;
463  size_t i, m=N%4;
464  for(i=0;i<m;i++) {
465  temp=v1[i];
466  v1[i]=v2[i];
467  v2[i]=temp;
468  }
469  for(i=m;i+3<N;i+=4) {
470  temp=v1[i];
471  v1[i]=v2[i];
472  v2[i]=temp;
473  temp=v1[i+1];
474  v1[i+1]=v2[i+1];
475  v2[i+1]=temp;
476  temp=v1[i+2];
477  v1[i+2]=v2[i+2];
478  v2[i+2]=temp;
479  temp=v1[i+3];
480  v1[i+3]=v2[i+3];
481  v2[i+3]=temp;
482  }
483  return;
484  }
485 
486  /** \brief Swap of of the first N elements of two
487  double-precision vectors
488 
489  This function swaps the elements of \c v1 and \c v2, one element
490  at a time.
491  */
492  template<class vec_t, class vec2_t>
493  void vector_swap_double(size_t N, vec_t &v1, vec2_t &v2) {
494  return vector_swap<vec_t,vec2_t,double>(N,v1,v2);
495  }
496 
497  /** \brief Swap of all the elements in two
498  double-precision vectors
499 
500  This function swaps the elements of \c v1 and \c v2, one element
501  at a time.
502 
503  \note It is almost always better to use <tt>std::swap</tt>
504  than this function, which is provided only in cases where
505  one knows one is going to be forced to use a vector type
506  without a properly defined <tt>std::swap</tt> method.
507  */
508  template<class vec_t, class vec2_t>
509  void vector_swap_double(vec_t &v1, vec2_t &v2) {
510  return vector_swap<vec_t,vec2_t,double>(v1,v2);
511  }
512 
513  /** \brief Swap two elements in a vector
514 
515  This function swaps the element \c i and element \c j of vector
516  \c v1.
517  */
518  template<class vec_t, class data_t>
519  void vector_swap(vec_t &v, size_t i, size_t j) {
520  data_t temp=v[i];
521  v[i]=v[j];
522  v[j]=temp;
523  return;
524  }
525 
526  /** \brief Swap two elements in a double-precision vector
527 
528  This function swaps the element \c i and element \c j of vector
529  \c v1.
530 
531  This function is used in \ref o2scl_linalg::QRPT_decomp() .
532  */
533  template<class vec_t>
534  void vector_swap_double(vec_t &v, size_t i, size_t j) {
535  return vector_swap<vec_t,double>(v,i,j);
536  }
537 
538  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
539  matrices
540 
541  This function swaps the elements of \c v1 and \c v2, one element
542  at a time.
543  */
544  template<class mat_t, class mat2_t, class data_t>
545  void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2) {
546  data_t temp;
547  for(size_t i=0;i<M;i++) {
548  for(size_t j=0;j<N;j++) {
549  temp=v1[i][j];
550  v1[i][j]=v2[i][j];
551  v2[i][j]=temp;
552  }
553  }
554  return;
555  }
556 
557  /** \brief Swap of the first \f$ (M,N) \f$ elements in two
558  double-precision matrices
559 
560  This function swaps the elements of \c m1 and \c m2, one element
561  at a time.
562  */
563  template<class mat_t, class mat2_t, class data_t>
564  void matrix_swap_double(size_t M, size_t N, mat_t &m1, mat2_t &m2) {
565  return matrix_swap<mat_t,mat2_t,double>(M,N,m1,m2);
566  }
567 
568  /** \brief Swap two elements in a matrix
569 
570  This function swaps the element <tt>(i1,j1)</tt> and
571  element <tt>(i2,j2)</tt> of matrix \c m1.
572  */
573  template<class mat_t, class data_t>
574  void matrix_swap(mat_t &m, size_t i1, size_t j1, size_t i2, size_t j2) {
575  data_t temp=m(i1,j1);
576  m(i1,j1)=m(i2,j2);
577  m(i2,j2)=temp;
578  return;
579  }
580 
581  /** \brief Swap two elements in a double-precision matrix
582 
583  This function swaps the element \c i and element \c j of matrix
584  \c v1.
585  */
586  template<class mat_t>
587  void matrix_swap_double(mat_t &m, size_t i1, size_t j1,
588  size_t i2, size_t j2) {
589  return matrix_swap<mat_t,double>(m,i1,j1,i2,j2);
590  }
591 
592  /** \brief Swap the first \c M rows of two columns in a matrix
593 
594  This function swaps the element <tt>(i1,j1)</tt> and
595  element <tt>(i2,j2)</tt> of matrix \c m1.
596  */
597  template<class mat_t, class data_t>
598  void matrix_swap_cols(size_t M, mat_t &m, size_t j1, size_t j2) {
599  data_t temp;
600  for(size_t i=0;i<M;i++) {
601  temp=m(i,j1);
602  m(i,j1)=m(i,j2);
603  m(i,j2)=temp;
604  }
605  return;
606  }
607 
608  /** \brief Swap the first \c M rows of two columns in a
609  double-precision matrix
610 
611  This function swaps the element \c i and element \c j of matrix
612  \c v1.
613  */
614  template<class mat_t>
615  void matrix_swap_cols_double(size_t M, mat_t &m, size_t j1, size_t j2) {
616  return matrix_swap_cols<mat_t,double>(M,m,j1,j2);
617  }
618 
619  /** \brief Swap the first \c N columns of two rows in a
620  matrix
621 
622  This function swaps the element <tt>(i1,j1)</tt> and
623  element <tt>(i2,j2)</tt> of matrix \c m1.
624  */
625  template<class mat_t, class data_t>
626  void matrix_swap_rows(size_t N, mat_t &m, size_t i1, size_t i2) {
627  data_t temp;
628  for(size_t j=0;j<N;j++) {
629  temp=m(i1,j);
630  m(i1,j)=m(i2,j);
631  m(i2,j)=temp;
632  }
633  return;
634  }
635 
636  /** \brief Swap the first \c N columns of two rows in a
637  double-precision matrix
638 
639  This function swaps the element \c i and element \c j of matrix
640  \c v1.
641  */
642  template<class mat_t>
643  void matrix_swap_rows_double(size_t N, mat_t &m, size_t i1, size_t i2) {
644  return matrix_swap_rows<mat_t,double>(N,m,i1,i2);
645  }
646  //@}
647 
648  /// \name Sorting vectors
649  //@{
650  /** \brief Provide a downheap() function for vector_sort()
651  */
652  template<class vec_t, class data_t>
653  void sort_downheap(vec_t &data, size_t n, size_t k) {
654 
655  data_t v=data[k];
656 
657  while (k<=n/2) {
658  size_t j=2*k;
659 
660  if (j<n && data[j] < data[j+1]) j++;
661  if (!(v < data[j])) break;
662  data[k]=data[j];
663  k=j;
664  }
665 
666  data[k]=v;
667 
668  return;
669  }
670 
671  /** \brief Sort a vector (in increasing order)
672 
673  This is a generic sorting template function using a heapsort
674  algorithm. It will work for any types \c data_t and \c vec_t for
675  which
676  - \c data_t has a non-const version of <tt>operator=</tt>
677  - \c data_t has a less than operator to compare elements
678  - <tt>vec_t::operator[]</tt> returns a non-const reference
679  to an object of type \c data_t
680 
681  In particular, it will work with the STL template class
682  <tt>std::vector</tt>, and arrays and pointers of numeric,
683  character, and string objects.
684 
685  For example,
686  \code
687  std::string list[3]={"dog","cat","fox"};
688  vector_sort<std::string[3],std::string>(3,list);
689  \endcode
690 
691  \note With this function template alone, the user cannot avoid
692  explicitly specifying the template types for this function
693  because there is no parameter of type \c data_t, and function
694  templates cannot handle default template types. For this
695  reason, the function template \ref o2scl::vector_sort_double() was
696  also created which provides the convenience of not requiring
697  the user to specify the vector template type.
698 
699  \note This sorting routine is not stable, i.e. equal elements
700  have arbtrary final ordering
701 
702  \note If \c n is zero, this function will do nothing and will
703  not call the error handler.
704 
705  This works similarly to the GSL function <tt>gsl_sort_vector()</tt>.
706  */
707  template<class vec_t, class data_t>
708  void vector_sort(size_t n, vec_t &data) {
709 
710  size_t N;
711  size_t k;
712 
713  if (n==0) return;
714 
715  N=n-1;
716  k=N/2;
717  k++;
718  do {
719  k--;
720  sort_downheap<vec_t,data_t>(data,N,k);
721  } while (k > 0);
722 
723  while (N > 0) {
724  data_t tmp=data[0];
725  data[0]=data[N];
726  data[N]=tmp;
727  N--;
728  sort_downheap<vec_t,data_t>(data,N,0);
729  }
730 
731  return;
732  }
733 
734  /** \brief Provide a downheap() function for vector_sort_index()
735  */
736  template<class vec_t, class vec_size_t>
737  void sort_index_downheap(size_t N, const vec_t &data, vec_size_t &order,
738  size_t k) {
739 
740  const size_t pki = order[k];
741 
742  while (k <= N / 2) {
743  size_t j = 2 * k;
744 
745  if (j < N && data[order[j]] < data[order[j + 1]]) {
746  j++;
747  }
748 
749  // [GSL] Avoid infinite loop if nan
750  if (!(data[pki] < data[order[j]])) {
751  break;
752  }
753 
754  order[k] = order[j];
755 
756  k = j;
757  }
758 
759  order[k] = pki;
760 
761  return;
762  }
763 
764  /** \brief Create a permutation which sorts a vector (in increasing order)
765 
766  This function takes a vector \c data and arranges a list of
767  indices in \c order, which give a sorted version of the vector.
768  The value <tt>order[i]</tt> gives the index of entry in in \c
769  data which corresponds to the <tt>i</tt>th value in the sorted
770  vector. The vector \c data is unchanged by this function, and
771  the initial values in \c order are ignored. Before calling this
772  function, \c order must already be allocated as a vector of size
773  \c n.
774 
775  For example, after calling this function, a sorted version the
776  vector can be output with
777  \code
778  size_t n=5;
779  double data[5]={3.1,4.1,5.9,2.6,3.5};
780  permutation order(n);
781  vector_sort_index(n,data,order);
782  for(size_t i=0;i<n;i++) {
783  cout << data[order[i]] << endl;
784  }
785  \endcode
786 
787  To create a permutation which stores as its <tt>i</tt>th element,
788  the index of <tt>data[i]</tt> in the sorted vector, you can
789  invert the permutation created by this function.
790 
791  This is a generic sorting template function. It will work for
792  any types \c vec_t and \c vec_size_t for which
793  - \c vec_t has an <tt>operator[]</tt>, and
794  - \c vec_size_t has an <tt>operator[]</tt> which returns
795  a \c size_t .
796  One possible type for \c vec_size_t is \ref o2scl::permutation.
797 
798  This works similarly to the GSL function <tt>gsl_sort_index()</tt>.
799  */
800  template<class vec_t, class vec_size_t>
801  void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order) {
802  size_t N;
803  size_t i, k;
804 
805  if (n == 0) return;
806 
807  // [GSL] Set permutation to identity
808 
809  for (i = 0 ; i < n ; i++) {
810  order[i] = i;
811  }
812 
813  /* [GSL] We have n_data elements, last element is at 'n_data-1',
814  first at '0' Set N to the last element number.
815  */
816  N = n - 1;
817 
818  k = N / 2;
819  // [GSL] Compensate the first use of 'k--'
820  k++;
821  do {
822  k--;
823  sort_index_downheap<vec_t,vec_size_t>(N,data,order,k);
824  } while (k > 0);
825 
826  while (N > 0) {
827 
828  // [GSL] First swap the elements
829  size_t tmp = order[0];
830  order[0] = order[N];
831  order[N] = tmp;
832 
833  // [GSL] Then process the heap
834  N--;
835 
836  sort_index_downheap<vec_t,vec_size_t>(N,data,order,0);
837  }
838 
839  return;
840  }
841 
842  /** \brief Sort a vector of doubles (in increasing order)
843 
844  This function is just a wrapper for
845  \code
846  vector_sort<vec_t,double>(n,data);
847  \endcode
848  See the documentation of \ref o2scl::vector_sort() for more
849  details.
850  */
851  template<class vec_t>
852  void vector_sort_double(size_t n, vec_t &data) {
853  return vector_sort<vec_t,double>(n,data);
854  }
855  //@}
856 
857  /// \name Smallest or largest subset functions
858  //@{
859  /** \brief Find the k smallest entries of the first \c n elements
860  of a vector
861 
862  Given a vector \c data of size \c n, this function sets the
863  first \c k entries of the vector \c smallest to the k smallest
864  entries from vector \c data in ascending order. The vector \c
865  smallest must be allocated beforehand to hold at least \c k
866  elements.
867 
868  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
869 
870  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
871  \f$ k << N \f$.
872 
873  If \c k is zero, then this function does nothing and
874  returns \ref o2scl::success .
875  */
876  template<class vec_t, class data_t>
877  void vector_smallest(size_t n, vec_t &data, size_t k, vec_t &smallest) {
878  if (k>n) {
879  O2SCL_ERR2("Subset length greater than size in ",
880  "function vector_smallest().",exc_einval);
881  }
882  if (k==0 || n==0) {
883  O2SCL_ERR2("Vector size zero or k zero in ",
884  "function vector_smallest().",exc_einval);
885  }
886 
887  // Take the first element
888  size_t j=1;
889  data_t xbound=data[0];
890  smallest[0]=xbound;
891 
892  // Examine the remaining elements
893  for(size_t i=1;i<n;i++) {
894  data_t xi=data[i];
895  if (j<k) {
896  j++;
897  } else if (xi>=xbound) {
898  continue;
899  }
900  size_t i1;
901  for(i1=j-1;i1>0;i1--) {
902  if (xi>smallest[i1-1]) break;
903  smallest[i1]=smallest[i1-1];
904  }
905  smallest[i1]=xi;
906  xbound=smallest[j-1];
907  }
908  return;
909  }
910 
911  /** \brief Find the k smallest entries of a vector
912  of a vector
913 
914  Given a vector \c data, this function sets the first \c k
915  entries of the vector \c smallest to the k smallest entries from
916  vector \c data in ascending order. The vector \c smallest
917  is resized if necessary to hold at least \c k elements.
918 
919  This works similarly to the GSL function <tt>gsl_sort_smallest()</tt>.
920 
921  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
922  \f$ k << N \f$.
923 
924  If \c k is zero, then this function does nothing and
925  returns \ref o2scl::success .
926  */
927  template<class vec_t, class data_t>
928  void vector_smallest(vec_t &data, size_t k, vec_t &smallest) {
929  size_t n=data.size();
930  if (smallest.size()<k) smallest.resize(k);
931  return vector_smallest(n,data,k,smallest);
932  }
933 
934  /** \brief Find the indexes of the k smallest entries among the
935  first \c n entries of a vector
936 
937  Given a vector \c data, this function sets the first \c k
938  entries of the vector \c smallest equal to the indexes of the k
939  smallest entries from vector \c data in ascending order. The
940  vector \c smallest is resized if necessary to hold at least \c k
941  elements.
942 
943  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
944  \f$ k << N \f$.
945 
946  If \c k is zero or \c n is zero or \f$ k > n\f$, then this
947  function calls the error handler.
948  */
949  template<class vec_t, class data_t, class vec_size_t>
950  void vector_smallest_index(size_t n, vec_t &data, size_t k,
951  vec_size_t &index) {
952  if (k>n) {
953  O2SCL_ERR2("Subset length greater than size in ",
954  "function vector_smallest_index().",exc_einval);
955  }
956  if (k==0 || n==0) {
957  O2SCL_ERR2("Vector size zero or k zero in ",
958  "function vector_smallest_index().",exc_einval);
959  }
960 
961  index.resize(k);
962 
963  // [GSL] Take the first element
964  size_t j=1;
965  data_t xbound=data[0];
966  index[0]=0;
967 
968  // [GSL] Examine the remaining elements
969  for(size_t i=1;i<n;i++) {
970  data_t xi=data[i];
971  if (j<k) {
972  j++;
973  } else if (xi>=xbound) {
974  continue;
975  }
976  size_t i1;
977  for(i1=j-1;i1>0;i1--) {
978  if (xi>data[index[i1-1]]) break;
979  index[i1]=index[i1-1];
980  }
981  index[i1]=i;
982  xbound=data[index[j-1]];
983  }
984  return;
985  }
986 
987  /** \brief Find the indexes of the k smallest entries of a vector
988  */
989  template<class vec_t, class data_t, class vec_size_t>
990  void vector_smallest_index(vec_t &data, size_t k,
991  vec_size_t &index) {
992  size_t n=data.size();
993  if (index.size()<k) index.resize(k);
994  return o2scl::vector_smallest_index<vec_t,data_t,vec_size_t>
995  (n,data,k,index);
996  }
997 
998  /** \brief Find the k largest entries of the first \c n elements
999  of a vector
1000 
1001  Given a vector \c data of size \c n this sets the first \c k
1002  entries of the vector \c largest to the k largest entries from
1003  vector \c data in descending order. The vector \c largest must
1004  be allocated beforehand to hold at least \c k elements.
1005 
1006  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1007 
1008  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1009  \f$ k << N \f$.
1010 
1011  If \c k is zero, then this function does nothing and
1012  returns \ref o2scl::success .
1013  */
1014  template<class vec_t, class data_t>
1015  void vector_largest(size_t n, vec_t &data, size_t k, vec_t &largest) {
1016  if (k>n) {
1017  O2SCL_ERR2("Subset length greater than size in ",
1018  "function vector_largest().",exc_einval);
1019  }
1020  if (k==0 || n==0) {
1021  O2SCL_ERR2("Vector size zero or k zero in ",
1022  "function vector_largest().",exc_einval);
1023  }
1024 
1025  // Take the first element
1026  size_t j=1;
1027  data_t xbound=data[0];
1028  largest[0]=xbound;
1029 
1030  // Examine the remaining elements
1031  for(size_t i=1;i<n;i++) {
1032  data_t xi=data[i];
1033  if (j<k) {
1034  j++;
1035  } else if (xi<=xbound) {
1036  continue;
1037  }
1038  size_t i1;
1039  for(i1=j-1;i1>0;i1--) {
1040  if (xi<largest[i1-1]) break;
1041  largest[i1]=largest[i1-1];
1042  }
1043  largest[i1]=xi;
1044  xbound=largest[j-1];
1045  }
1046  return;
1047  }
1048 
1049  /** \brief Find the k largest entries of a vector
1050  of a vector
1051 
1052  Given a vector \c data, this function sets the first \c k
1053  entries of the vector \c largest to the k largest entries from
1054  vector \c data in ascending order. The vector \c largest
1055  is resized if necessary to hold at least \c k elements.
1056 
1057  This works similarly to the GSL function <tt>gsl_sort_largest()</tt>.
1058 
1059  \note This \f$ {\cal O}(k N) \f$ algorithm is useful only when
1060  \f$ k << N \f$.
1061 
1062  If \c k is zero, then this function does nothing and
1063  returns \ref o2scl::success .
1064  */
1065  template<class vec_t, class data_t>
1066  void vector_largest(vec_t &data, size_t k, vec_t &largest) {
1067  size_t n=data.size();
1068  if (largest.size()<k) largest.resize(k);
1069  return vector_largest(n,data,k,largest);
1070  }
1071  //@}
1072 
1073  /// \name Vector minimum and maximum functions
1074  //@{
1075  /** \brief Compute the maximum of the first \c n elements of a vector
1076  */
1077  template<class vec_t, class data_t>
1078  data_t vector_max_value(size_t n, const vec_t &data) {
1079 
1080  if (n==0) {
1081  O2SCL_ERR("Sent size=0 to vector_max_value().",exc_efailed);
1082  }
1083  data_t max=data[0];
1084  for(size_t i=1;i<n;i++) {
1085  if (data[i]>max) {
1086  max=data[i];
1087  }
1088  }
1089  return max;
1090  }
1091 
1092  /** \brief Compute the maximum value of a vector
1093  */
1094  template<class vec_t, class data_t>
1095  data_t vector_max_value(const vec_t &data) {
1096 
1097  size_t n=data.size();
1098  if (n==0) {
1099  O2SCL_ERR("Sent empty vector to vector_max_value().",exc_efailed);
1100  }
1101  data_t max=data[0];
1102  for(size_t i=1;i<n;i++) {
1103  if (data[i]>max) {
1104  max=data[i];
1105  }
1106  }
1107  return max;
1108  }
1109 
1110  /** \brief Compute the index which holds the
1111  maximum of the first \c n elements of a vector
1112  */
1113  template<class vec_t, class data_t>
1114  size_t vector_max_index(size_t n, const vec_t &data) {
1115 
1116  if (n==0) {
1117  O2SCL_ERR("Sent size=0 to vector_max_index().",exc_efailed);
1118  }
1119  data_t max=data[0];
1120  size_t ix=0;
1121  for(size_t i=1;i<n;i++) {
1122  if (data[i]>max) {
1123  max=data[i];
1124  ix=i;
1125  }
1126  }
1127  return ix;
1128  }
1129 
1130  /** \brief Compute the maximum of the first \c n elements of a vector
1131  */
1132  template<class vec_t, class data_t>
1133  void vector_max(size_t n, const vec_t &data, size_t &index,
1134  data_t &val) {
1135 
1136  if (n==0) {
1137  O2SCL_ERR("Sent size=0 to vector_max().",exc_efailed);
1138  }
1139  val=data[0];
1140  index=0;
1141  for(size_t i=1;i<n;i++) {
1142  if (data[i]>val) {
1143  val=data[i];
1144  index=i;
1145  }
1146  }
1147  return;
1148  }
1149 
1150  /** \brief Compute the minimum of the first \c n elements of a vector
1151  */
1152  template<class vec_t, class data_t>
1153  data_t vector_min_value(size_t n, const vec_t &data) {
1154 
1155  if (n==0) {
1156  O2SCL_ERR("Sent size=0 to vector_min_value().",exc_efailed);
1157  }
1158  data_t min=data[0];
1159  for(size_t i=1;i<n;i++) {
1160  if (data[i]<min) {
1161  min=data[i];
1162  }
1163  }
1164  return min;
1165  }
1166 
1167  /** \brief Compute the minimum value in a vector
1168  */
1169  template<class vec_t, class data_t>
1170  data_t vector_min_value(const vec_t &data) {
1171 
1172  size_t n=data.size();
1173  if (n==0) {
1174  O2SCL_ERR("Sent empty vector to vector_min_value().",exc_efailed);
1175  }
1176  data_t min=data[0];
1177  for(size_t i=1;i<n;i++) {
1178  if (data[i]<min) {
1179  min=data[i];
1180  }
1181  }
1182  return min;
1183  }
1184 
1185  /** \brief Compute the index which holds the
1186  minimum of the first \c n elements of a vector
1187  */
1188  template<class vec_t, class data_t>
1189  size_t vector_min_index(size_t n, const vec_t &data) {
1190 
1191  if (n==0) {
1192  O2SCL_ERR("Sent size=0 to vector_min_index().",exc_efailed);
1193  }
1194  data_t min=data[0];
1195  size_t ix=0;
1196  for(size_t i=1;i<n;i++) {
1197  if (data[i]<min) {
1198  min=data[i];
1199  ix=i;
1200  }
1201  }
1202  return ix;
1203  }
1204 
1205  /** \brief Compute the minimum of the first \c n elements of a vector
1206  */
1207  template<class vec_t, class data_t>
1208  void vector_min(size_t n, const vec_t &data, size_t &index,
1209  data_t &val) {
1210 
1211  if (n==0) {
1212  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1213  }
1214  val=data[0];
1215  index=0;
1216  for(size_t i=1;i<n;i++) {
1217  if (data[i]<val) {
1218  val=data[i];
1219  index=i;
1220  }
1221  }
1222  return;
1223  }
1224 
1225  /** \brief Compute the minimum and maximum of the first
1226  \c n elements of a vector
1227  */
1228  template<class vec_t, class data_t>
1229  void vector_minmax_value(size_t n, vec_t &data,
1230  data_t &min, data_t &max) {
1231 
1232  if (n==0) {
1233  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1234  }
1235  min=data[0];
1236  max=min;
1237  for(size_t i=1;i<n;i++) {
1238  if (data[i]<min) {
1239  min=data[i];
1240  }
1241  if (data[i]>max) {
1242  max=data[i];
1243  }
1244  }
1245  return;
1246  }
1247 
1248  /** \brief Compute the minimum and maximum of the first
1249  \c n elements of a vector
1250  */
1251  template<class vec_t, class data_t>
1252  void vector_minmax_index(size_t n, vec_t &data,
1253  size_t &ix_min, size_t &ix_max) {
1254 
1255  if (n==0) {
1256  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1257  }
1258  data_t min=data[0];
1259  data_t max=min;
1260  ix_min=0;
1261  ix_max=0;
1262  for(size_t i=1;i<n;i++) {
1263  if (data[i]<min) {
1264  min=data[i];
1265  ix_min=i;
1266  }
1267  if (data[i]>max) {
1268  max=data[i];
1269  ix_max=i;
1270  }
1271  }
1272  return;
1273  }
1274 
1275  /** \brief Compute the minimum and maximum of the first
1276  \c n elements of a vector
1277  */
1278  template<class vec_t, class data_t>
1279  void vector_minmax(size_t n, vec_t &data,
1280  size_t &ix_min, data_t &min,
1281  size_t &ix_max, data_t &max) {
1282 
1283  if (n==0) {
1284  O2SCL_ERR("Sent size=0 to vector_min().",exc_efailed);
1285  }
1286  min=data[0];
1287  max=min;
1288  ix_min=0;
1289  ix_max=0;
1290  for(size_t i=1;i<n;i++) {
1291  if (data[i]<min) {
1292  min=data[i];
1293  ix_min=i;
1294  }
1295  if (data[i]>max) {
1296  max=data[i];
1297  ix_max=i;
1298  }
1299  }
1300  return;
1301  }
1302  //@}
1303 
1304  /// \name Minima and maxima of vectors through quadratic fit
1305  //@{
1306  /** \brief Maximum of vector by quadratic fit
1307  */
1308  template<class vec_t, class data_t>
1309  data_t vector_max_quad(size_t n, const vec_t &data) {
1310  size_t ix=vector_max_index<vec_t,data_t>(n,data);
1311  if (ix==0) {
1312  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1313  } else if (ix==n-1) {
1314  return quadratic_extremum_y<data_t>
1315  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1316  }
1317  return quadratic_extremum_y<data_t>
1318  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1319  }
1320 
1321  /** \brief Maximum of vector by quadratic fit
1322  */
1323  template<class vec_t, class data_t>
1324  data_t vector_max_quad(size_t n, const vec_t &x, const vec_t &y) {
1325  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1326  if (ix==0 || ix==n-1) return y[ix];
1327  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1328  y[ix-1],y[ix],y[ix+1]);
1329  }
1330 
1331  /** \brief Location of vector maximum by quadratic fit
1332  */
1333  template<class vec_t, class data_t>
1334  data_t vector_max_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1335  size_t ix=vector_max_index<vec_t,data_t>(n,y);
1336  if (ix==0 || ix==n-1) return y[ix];
1337  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1338  y[ix-1],y[ix],y[ix+1]);
1339  }
1340 
1341  /** \brief Minimum of vector by quadratic fit
1342  */
1343  template<class vec_t, class data_t>
1344  data_t vector_min_quad(size_t n, const vec_t &data) {
1345  size_t ix=vector_min_index<vec_t,data_t>(n,data);
1346  if (ix==0) {
1347  return quadratic_extremum_y<data_t>(0,1,2,data[0],data[1],data[2]);
1348  } else if (ix==n-1) {
1349  return quadratic_extremum_y<data_t>
1350  (n-3,n-2,n-1,data[n-3],data[n-2],data[n-1]);
1351  }
1352  return quadratic_extremum_y<data_t>
1353  (ix-1,ix,ix+1,data[ix-1],data[ix],data[ix+1]);
1354  }
1355 
1356  /** \brief Minimum of vector by quadratic fit
1357  */
1358  template<class vec_t, class data_t>
1359  data_t vector_min_quad(size_t n, const vec_t &x, const vec_t &y) {
1360  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1361  if (ix==0 || ix==n-1) return y[ix];
1362  return quadratic_extremum_y<data_t>(x[ix-1],x[ix],x[ix+1],
1363  y[ix-1],y[ix],y[ix+1]);
1364  }
1365 
1366  /** \brief Location of vector minimum by quadratic fit
1367  */
1368  template<class vec_t, class data_t>
1369  data_t vector_min_quad_loc(size_t n, const vec_t &x, const vec_t &y) {
1370  size_t ix=vector_min_index<vec_t,data_t>(n,y);
1371  if (ix==0 || ix==n-1) return y[ix];
1372  return quadratic_extremum_x<data_t>(x[ix-1],x[ix],x[ix+1],
1373  y[ix-1],y[ix],y[ix+1]);
1374  }
1375  //@}
1376 
1377  /// \name Matrix minimum and maximum functions
1378  //@{
1379  /** \brief Compute the maximum of the lower-left part of a matrix
1380  */
1381  template<class mat_t, class data_t>
1382  data_t matrix_max_value(size_t m, const size_t n, const mat_t &data) {
1383 
1384  if (m==0 || n==0) {
1385  std::string str=((std::string)"Matrix with zero size (")+
1386  o2scl::itos(m)+","+o2scl::itos(n)+") in "+
1387  "matrix_max_value().";
1388  O2SCL_ERR(str.c_str(),exc_einval);
1389  }
1390  data_t max=data(0,0);
1391  for(size_t i=0;i<m;i++) {
1392  for(size_t j=0;j<n;j++) {
1393  if (data(i,j)>max) {
1394  max=data(i,j);
1395  }
1396  }
1397  }
1398  return max;
1399  }
1400 
1401  /** \brief Compute the maximum of a matrix
1402  */
1403  template<class mat_t, class data_t> data_t
1404  matrix_max_value(const mat_t &data) {
1405  size_t m=data.size1();
1406  size_t n=data.size2();
1407  if (m==0 || n==0) {
1408  std::string str=((std::string)"Matrix with zero size (")+
1409  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1410  "matrix_max_value().";
1411  O2SCL_ERR(str.c_str(),exc_einval);
1412  }
1413  data_t max=data(0,0);
1414  for(size_t i=0;i<m;i++) {
1415  for(size_t j=0;j<n;j++) {
1416  if (data(i,j)>max) {
1417  max=data(i,j);
1418  }
1419  }
1420  }
1421  return max;
1422  }
1423 
1424  /** \brief Compute the maximum of a matrix
1425  */
1426  template<class mat_t> double
1427  matrix_max_value_double(const mat_t &data) {
1428  size_t m=data.size1();
1429  size_t n=data.size2();
1430  if (m==0 || n==0) {
1431  std::string str=((std::string)"Matrix with zero size (")+
1432  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1433  "matrix_max_value_double().";
1434  O2SCL_ERR(str.c_str(),exc_einval);
1435  }
1436  double max=data(0,0);
1437  for(size_t i=0;i<m;i++) {
1438  for(size_t j=0;j<n;j++) {
1439  if (data(i,j)>max) {
1440  max=data(i,j);
1441  }
1442  }
1443  }
1444  return max;
1445  }
1446 
1447  /** \brief Compute the maximum of a matrix and return
1448  the indices of the maximum element
1449  */
1450  template<class mat_t, class data_t>
1451  void matrix_max_index(size_t m, size_t n, const mat_t &data,
1452  size_t &i_max, size_t &j_max, data_t &max) {
1453 
1454  if (m==0 || n==0) {
1455  std::string str=((std::string)"Matrix with zero size (")+
1456  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1457  "matrix_max_index().";
1458  O2SCL_ERR(str.c_str(),exc_einval);
1459  }
1460  max=data(0,0);
1461  i_max=0;
1462  j_max=0;
1463  for(size_t i=0;i<m;i++) {
1464  for(size_t j=0;j<n;j++) {
1465  if (data(i,j)>max) {
1466  max=data(i,j);
1467  i_max=i;
1468  j_max=j;
1469  }
1470  }
1471  }
1472  return;
1473  }
1474 
1475  /** \brief Compute the maximum of a matrix and return
1476  the indices of the maximum element
1477  */
1478  template<class mat_t, class data_t>
1479  void matrix_max_index(const mat_t &data,
1480  size_t &i_max, size_t &j_max, data_t &max) {
1481 
1482  size_t m=data.size1();
1483  size_t n=data.size2();
1484  if (m==0 || n==0) {
1485  std::string str=((std::string)"Matrix with zero size (")+
1486  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1487  "matrix_max_index().";
1488  O2SCL_ERR(str.c_str(),exc_einval);
1489  }
1490  max=data(0,0);
1491  i_max=0;
1492  j_max=0;
1493  for(size_t i=0;i<m;i++) {
1494  for(size_t j=0;j<n;j++) {
1495  if (data(i,j)>max) {
1496  max=data(i,j);
1497  i_max=i;
1498  j_max=j;
1499  }
1500  }
1501  }
1502  return;
1503  }
1504 
1505  /** \brief Compute the minimum of a matrix
1506  */
1507  template<class mat_t, class data_t>
1508  data_t matrix_min_value(size_t m, size_t n, const mat_t &data) {
1509 
1510  if (m==0 || n==0) {
1511  std::string str=((std::string)"Matrix with zero size (")+
1512  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1513  "matrix_min_value().";
1514  O2SCL_ERR(str.c_str(),exc_einval);
1515  }
1516  data_t min=data(0,0);
1517  for(size_t i=0;i<m;i++) {
1518  for(size_t j=0;j<n;j++) {
1519  if (data(i,j)<min) {
1520  min=data(i,j);
1521  }
1522  }
1523  }
1524  return min;
1525  }
1526 
1527  /** \brief Compute the minimum of a matrix
1528  */
1529  template<class mat_t, class data_t>
1530  data_t matrix_min_value(const mat_t &data) {
1531 
1532  size_t m=data.size1();
1533  size_t n=data.size2();
1534  if (m==0 || n==0) {
1535  std::string str=((std::string)"Matrix with zero size (")+
1536  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1537  "matrix_min_value().";
1538  O2SCL_ERR(str.c_str(),exc_einval);
1539  }
1540  data_t min=data(0,0);
1541  for(size_t i=0;i<m;i++) {
1542  for(size_t j=0;j<n;j++) {
1543  if (data(i,j)<min) {
1544  min=data(i,j);
1545  }
1546  }
1547  }
1548  return min;
1549  }
1550 
1551  /** \brief Compute the minimum of a matrix
1552  */
1553  template<class mat_t>
1554  double matrix_min_value_double(const mat_t &data) {
1555 
1556  size_t m=data.size1();
1557  size_t n=data.size2();
1558  if (m==0 || n==0) {
1559  std::string str=((std::string)"Matrix with zero size (")+
1560  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1561  "matrix_min_value().";
1562  O2SCL_ERR(str.c_str(),exc_einval);
1563  }
1564  double min=data(0,0);
1565  for(size_t i=0;i<m;i++) {
1566  for(size_t j=0;j<n;j++) {
1567  if (data(i,j)<min) {
1568  min=data(i,j);
1569  }
1570  }
1571  }
1572  return min;
1573  }
1574 
1575  /** \brief Compute the minimum of a matrix and return
1576  the indices of the minimum element
1577  */
1578  template<class mat_t, class data_t>
1579  void matrix_min_index(size_t n, size_t m, const mat_t &data,
1580  size_t &i_min, size_t &j_min, data_t &min) {
1581 
1582  if (m==0 || n==0) {
1583  std::string str=((std::string)"Matrix with zero size (")+
1584  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1585  "matrix_min_index().";
1586  O2SCL_ERR(str.c_str(),exc_einval);
1587  }
1588  min=data(0,0);
1589  i_min=0;
1590  j_min=0;
1591  for(size_t i=0;i<m;i++) {
1592  for(size_t j=0;j<n;j++) {
1593  if (data(i,j)<min) {
1594  min=data(i,j);
1595  i_min=i;
1596  j_min=j;
1597  }
1598  }
1599  }
1600  return;
1601  }
1602 
1603  /** \brief Compute the minimum of a matrix and return
1604  the indices of the minimum element
1605  */
1606  template<class mat_t, class data_t>
1607  void matrix_min_index(const mat_t &data,
1608  size_t &i_min, size_t &j_min, data_t &min) {
1609 
1610  size_t m=data.size1();
1611  size_t n=data.size2();
1612  if (m==0 || n==0) {
1613  std::string str=((std::string)"Matrix with zero size (")+
1614  o2scl::szttos(m)+","+o2scl::szttos(n)+") in "+
1615  "matrix_min_index().";
1616  O2SCL_ERR(str.c_str(),exc_einval);
1617  }
1618  min=data(0,0);
1619  i_min=0;
1620  j_min=0;
1621  for(size_t i=0;i<m;i++) {
1622  for(size_t j=0;j<n;j++) {
1623  if (data(i,j)<min) {
1624  min=data(i,j);
1625  i_min=i;
1626  j_min=j;
1627  }
1628  }
1629  }
1630  return;
1631  }
1632 
1633  /** \brief Compute the minimum and maximum of a matrix
1634  */
1635  template<class mat_t, class data_t>
1636  void matrix_minmax(size_t n, size_t m, const mat_t &data,
1637  data_t &min, data_t &max) {
1638 
1639  if (n==0 || m==0) {
1640  O2SCL_ERR("Sent size=0 to matrix_min().",exc_efailed);
1641  }
1642  min=data(0,0);
1643  max=data(0,0);
1644  for(size_t i=0;i<n;i++) {
1645  for(size_t j=0;j<m;j++) {
1646  if (data(i,j)<min) {
1647  min=data(i,j);
1648  } else if (data(i,j)>max) {
1649  max=data(i,j);
1650  }
1651  }
1652  }
1653  return;
1654  }
1655 
1656  /** \brief Compute the minimum and maximum of a matrix
1657  */
1658  template<class mat_t, class data_t>
1659  void matrix_minmax(const mat_t &data,
1660  data_t &min, data_t &max) {
1661  return matrix_minmax<mat_t,data_t>
1662  (data.size1(),data.size2(),data,min,max);
1663  }
1664 
1665  /** \brief Compute the minimum and maximum of a matrix and
1666  return their locations
1667  */
1668  template<class mat_t, class data_t>
1669  void matrix_minmax_index(size_t n, size_t m, const mat_t &data,
1670  size_t &i_min, size_t &j_min, data_t &min,
1671  size_t &i_max, size_t &j_max, data_t &max) {
1672 
1673  if (n==0 || m==0) {
1674  O2SCL_ERR2("Sent size=0 to function ",
1675  "matrix_minmax_index().",exc_efailed);
1676  }
1677  min=data(0,0);
1678  i_min=0;
1679  j_min=0;
1680  max=data(0,0);
1681  i_max=0;
1682  j_max=0;
1683  for(size_t i=0;i<n;i++) {
1684  for(size_t j=0;j<m;j++) {
1685  if (data(i,j)<min) {
1686  min=data(i,j);
1687  i_min=i;
1688  j_min=j;
1689  } else if (data(i,j)>max) {
1690  max=data(i,j);
1691  i_max=i;
1692  j_max=j;
1693  }
1694  }
1695  }
1696  return;
1697  }
1698 
1699  /** \brief Compute the sum of matrix elements
1700  */
1701  template<class mat_t, class data_t>
1702  data_t matrix_sum(size_t m, size_t n, const mat_t &data) {
1703 
1704  data_t sum=0.0;
1705  for(size_t i=0;i<m;i++) {
1706  for(size_t j=0;j<n;j++) {
1707  sum+=data(i,j);
1708  }
1709  }
1710  return sum;
1711  }
1712 
1713  /** \brief Compute the sum of matrix elements
1714  */
1715  template<class mat_t, class data_t>
1716  data_t matrix_sum(const mat_t &data) {
1717  return matrix_sum<mat_t,data_t>(data.size1(),data.size2(),data);
1718  }
1719  //@}
1720 
1721  /// \name Searching vectors and matrices
1722  //@{
1723  /** \brief Lookup the value \c x0 in the first \c n elements of
1724  vector \c x
1725 
1726  The function finds the element among the first \c n elements of
1727  \c x which is closest to the value \c x0. It ignores all
1728  elements in \c x which are not finite. If the vector is empty,
1729  or if all of the first \c n elements in \c x are not finite,
1730  then the error handler will be called.
1731 
1732  This function works for all classes \c vec_t where an operator[]
1733  is defined which returns a double (either as a value or a
1734  reference).
1735  */
1736  template<class vec_t>
1737  size_t vector_lookup(size_t n, const vec_t &x, double x0) {
1738  if (n==0) {
1739  O2SCL_ERR("Empty vector in function vector_lookup().",
1740  exc_einval);
1741  return 0;
1742  }
1743  size_t row=0, i=0;
1744  while(!std::isfinite(x[i]) && i<n-1) i++;
1745  if (i==n-1) {
1746  O2SCL_ERR2("Entire vector not finite in ",
1747  "function vector_lookup()",exc_einval);
1748  return 0;
1749  }
1750  double best=x[i], bdiff=fabs(x[i]-x0);
1751  for(;i<n;i++) {
1752  if (std::isfinite(x[i]) && fabs(x[i]-x0)<bdiff) {
1753  row=i;
1754  best=x[i];
1755  bdiff=fabs(x[i]-x0);
1756  }
1757  }
1758  return row;
1759  }
1760 
1761  /** \brief Lookup element \c x0 in vector \c x
1762 
1763  This function finds the element in vector \c x which is closest
1764  to \c x0. It ignores all elements in \c x which are not finite.
1765  If the vector is empty, or if all of the
1766  elements in \c x are not finite, then the error handler will be
1767  called.
1768 
1769  This function works for all classes \c vec_t with a
1770  <tt>size()</tt> method and where an operator[] is defined which
1771  returns a double (either as a value or a reference).
1772  */
1773  template<class vec_t>
1774  size_t vector_lookup(const vec_t &x, double x0) {
1775  return vector_lookup(x.size(),x,x0);
1776  }
1777 
1778  /** \brief Lookup an element in the first $(m,n)$ entries in a matrix
1779 
1780  Return the location <tt>(i,j)</tt> of the element closest to
1781  \c x0.
1782  */
1783  template<class mat_t>
1784  void matrix_lookup(size_t m, size_t n, const mat_t &A,
1785  double x0, size_t &i, size_t &j) {
1786  if (m==0 || n==0) {
1787  O2SCL_ERR("Empty matrix in matrix_lookup().",
1788  exc_einval);
1789  }
1790  double dist=0.0;
1791  bool found_one=false;
1792  for(size_t i2=0;i2<m;i2++) {
1793  for(size_t j2=0;j2<n;j2++) {
1794  if (std::isfinite(A(i,j))) {
1795  if (found_one==false) {
1796  dist=fabs(A(i,j)-x0);
1797  found_one=true;
1798  i=i2;
1799  j=j2;
1800  } else {
1801  if (fabs(A(i,j)-x0)<dist) {
1802  dist=fabs(A(i,j)-x0);
1803  i=i2;
1804  j=j2;
1805  }
1806  }
1807  }
1808  }
1809  }
1810  if (found_one==false) {
1811  O2SCL_ERR2("Entire matrix not finite in ",
1812  "function matrix_lookup()",exc_einval);
1813  }
1814  return;
1815  }
1816 
1817  /** \brief Lookup an element in a matrix
1818 
1819  Return the location <tt>(i,j)</tt> of the element closest to
1820  \c x0.
1821  */
1822  template<class mat_t>
1823  void matrix_lookup(const mat_t &A,
1824  double x0, size_t &i, size_t &j) {
1825  matrix_lookup(A.size1(),A.size2(),x0,i,j);
1826  return;
1827  }
1828 
1829  /** \brief Binary search a part of an increasing vector for
1830  <tt>x0</tt>.
1831 
1832  This function performs a binary search of between
1833  <tt>x[lo]</tt> and <tt>x[hi]</tt>. It returns
1834  - \c lo if \c x0 < <tt>x[lo]</tt>
1835  - \c i if <tt>x[i]</tt> <= \c x0 < <tt>x[i+2]</tt>
1836  for \c lo <= \c i < \c hi
1837  - \c hi-1 if \c x0 >= \c <tt>x[hi-1]</tt>
1838 
1839  This function is designed to find the interval containing \c x0,
1840  not the index of the element closest to x0. To perform the
1841  latter operation, you can use \ref vector_lookup().
1842 
1843  The element at <tt>x[hi]</tt> is never referenced by this
1844  function. The parameter \c hi can be either the index of the
1845  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1846  with starting index <tt>0</tt>), or the index of one element
1847  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1848  index <tt>0</tt>) for a depending on whether or not the user
1849  wants to allow the function to return the index of the last
1850  element.
1851 
1852  This function operates in the same way as
1853  <tt>gsl_interp_bsearch()</tt>.
1854 
1855  The operation of this function is undefined if the data is
1856  not strictly monotonic, i.e. if some of the data elements are
1857  equal.
1858 
1859  This function will call the error handler if \c lo is
1860  greater than \c hi.
1861  */
1862  template<class vec_t, class data_t>
1863  size_t vector_bsearch_inc(const data_t x0, const vec_t &x,
1864  size_t lo, size_t hi) {
1865  if (lo>hi) {
1866  O2SCL_ERR2("Low and high indexes backwards in ",
1867  "function vector_bsearch_inc().",exc_einval);
1868  }
1869  while (hi>lo+1) {
1870  size_t i=(hi+lo)/2;
1871  if (x[i]>x0) {
1872  hi=i;
1873  } else {
1874  lo=i;
1875  }
1876  }
1877 
1878  return lo;
1879  }
1880 
1881  /** \brief Binary search a part of an decreasing vector for
1882  <tt>x0</tt>.
1883 
1884  This function performs a binary search of between
1885  <tt>x[lo]</tt> and <tt>x[hi]</tt> (inclusive). It returns
1886  - \c lo if \c x0 > <tt>x[lo]</tt>
1887  - \c i if <tt>x[i]</tt> >= \c x0 > <tt>x[i+1]</tt>
1888  for \c lo <= \c i < \c hi
1889  - \c hi-1 if \c x0 <= \c <tt>x[hi-1]</tt>
1890 
1891  The element at <tt>x[hi]</tt> is never referenced by this
1892  function. The parameter \c hi can be either the index of the
1893  last element (e.g. <tt>n-1</tt> for a vector of size <tt>n</tt>
1894  with starting index <tt>0</tt>), or the index of one element
1895  (e.g. <tt>n</tt> for a vector of size <tt>n</tt> and a starting
1896  index <tt>0</tt>) for a depending on whether or not the user
1897  wants to allow the function to return the index of the last
1898  element.
1899 
1900  The operation of this function is undefined if the data is
1901  not strictly monotonic, i.e. if some of the data elements are
1902  equal.
1903 
1904  This function will call the error handler if \c lo is
1905  greater than \c hi.
1906  */
1907  template<class vec_t, class data_t>
1908  size_t vector_bsearch_dec(const data_t x0, const vec_t &x,
1909  size_t lo, size_t hi) {
1910  if (lo>hi) {
1911  O2SCL_ERR2("Low and high indexes backwards in ",
1912  "function vector_bsearch_dec().",exc_einval);
1913  }
1914  while (hi>lo+1) {
1915  size_t i=(hi+lo)/2;
1916  if (x[i]<x0) {
1917  hi=i;
1918  } else {
1919  lo=i;
1920  }
1921  }
1922 
1923  return lo;
1924  }
1925 
1926  /** \brief Binary search a part of a monotonic vector for
1927  <tt>x0</tt>.
1928 
1929  This wrapper just calls \ref o2scl::vector_bsearch_inc() or
1930  \ref o2scl::vector_bsearch_dec() depending on the ordering
1931  of \c x.
1932  */
1933  template<class vec_t, class data_t>
1934  size_t vector_bsearch(const data_t x0, const vec_t &x,
1935  size_t lo, size_t hi) {
1936  if (x[lo]<x[hi-1]) {
1937  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1938  }
1939  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
1940  }
1941 
1942  /** \brief Binary search a monotonic vector for
1943  <tt>x0</tt>.
1944 
1945  This function calls \ref o2scl::vector_bsearch_inc() or
1946  \ref o2scl::vector_bsearch_dec() depending on the ordering
1947  of \c x.
1948  */
1949  template<class vec_t, class data_t>
1950  size_t vector_bsearch(const data_t x0, const vec_t &x) {
1951  size_t lo=0;
1952  size_t hi=x.size();
1953  if (x[lo]<x[hi-1]) {
1954  return vector_bsearch_inc<vec_t,data_t>(x0,x,lo,hi);
1955  }
1956  return vector_bsearch_dec<vec_t,data_t>(x0,x,lo,hi);
1957  }
1958  //@}
1959 
1960  /// \name Miscellaneous mathematical functions
1961  //@{
1962  /** \brief Compute the sum of the first \c n elements of a vector
1963 
1964  If \c n is zero, this will return 0 without throwing
1965  an exception.
1966  */
1967  template<class vec_t, class data_t>
1968  data_t vector_sum(size_t n, vec_t &data) {
1969 
1970  data_t sum=0.0;
1971  for(size_t i=0;i<n;i++) {
1972  sum+=data[i];
1973  }
1974  return sum;
1975  }
1976 
1977  /** \brief Compute the sum of all the elements of a vector
1978 
1979  If the vector has zero size, this will return 0 without
1980  throwing an exception.
1981  */
1982  template<class vec_t, class data_t> data_t vector_sum(vec_t &data) {
1983  data_t sum=0.0;
1984  for(size_t i=0;i<data.size();i++) {
1985  sum+=data[i];
1986  }
1987  return sum;
1988  }
1989 
1990  /** \brief Compute the sum of the first \c n elements of a vector
1991  of double-precision numbers
1992 
1993  If \c n is zero, this will return 0 without throwing
1994  an exception.
1995  */
1996  template<class vec_t>double vector_sum_double(size_t n, vec_t &data) {
1997  double sum=0.0;
1998  for(size_t i=0;i<n;i++) {
1999  sum+=data[i];
2000  }
2001  return sum;
2002  }
2003 
2004  /** \brief Compute the sum of all the elements of a vector
2005  of double-precision numbers
2006 
2007  If the vector has zero size, this will return 0 without
2008  throwing an exception.
2009  */
2010  template<class vec_t> double vector_sum_double(vec_t &data) {
2011  double sum=0.0;
2012  for(size_t i=0;i<data.size();i++) {
2013  sum+=data[i];
2014  }
2015  return sum;
2016  }
2017 
2018  /** \brief Compute the norm of the first \c n entries of a
2019  vector of floating-point (single or double precision) numbers
2020 
2021  This function is a more generic version of
2022  \ref o2scl_cblas::dnrm2 .
2023  */
2024  template<class vec_t, class data_t>
2025  data_t vector_norm(size_t n, const vec_t &x) {
2026 
2027  data_t scale = 0.0;
2028  data_t ssq = 1.0;
2029 
2030  if (n <= 0) {
2031  return 0.0;
2032  } else if (n == 1) {
2033  return fabs(x[0]);
2034  }
2035 
2036  for (size_t i = 0; i < n; i++) {
2037  const data_t xx = x[i];
2038 
2039  if (xx != 0.0) {
2040  const data_t ax = fabs(xx);
2041 
2042  if (scale < ax) {
2043  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2044  scale = ax;
2045  } else {
2046  ssq += (ax / scale) * (ax / scale);
2047  }
2048  }
2049 
2050  }
2051 
2052  return scale * sqrt(ssq);
2053  }
2054 
2055  /** \brief Compute the norm of a vector of floating-point
2056  (single or double precision) numbers
2057  */
2058  template<class vec_t, class data_t> data_t vector_norm(const vec_t &x) {
2059  return vector_norm<vec_t,data_t>(x.size(),x);
2060  }
2061 
2062  /** \brief Compute the norm of the first \c n entries of a
2063  vector of double precision numbers
2064 
2065  This function is a more generic version of
2066  \ref o2scl_cblas::dnrm2 .
2067  */
2068  template<class vec_t>
2069  double vector_norm_double(size_t n, const vec_t &x) {
2070 
2071  double scale = 0.0;
2072  double ssq = 1.0;
2073 
2074  if (n <= 0) {
2075  return 0.0;
2076  } else if (n == 1) {
2077  return fabs(x[0]);
2078  }
2079 
2080  for (size_t i = 0; i < n; i++) {
2081  const double xx = x[i];
2082 
2083  if (xx != 0.0) {
2084  const double ax = fabs(xx);
2085 
2086  if (scale < ax) {
2087  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2088  scale = ax;
2089  } else {
2090  ssq += (ax / scale) * (ax / scale);
2091  }
2092  }
2093 
2094  }
2095 
2096  return scale * sqrt(ssq);
2097  }
2098 
2099  /** \brief Compute the norm of a vector of double precision numbers
2100  */
2101  template<class vec_t> double vector_norm_double(const vec_t &x) {
2102  return vector_norm_double<vec_t>(x.size(),x);
2103  }
2104  //@}
2105 
2106  /// \name Other vector and matrix functions
2107  //@{
2108  /** \brief "Rotate" a vector so that the kth element is now the beginning
2109 
2110  This is a generic template function which will work for
2111  any types \c data_t and \c vec_t for which
2112  - \c data_t has an <tt>operator=</tt>
2113  - <tt>vec_t::operator[]</tt> returns a reference
2114  to an object of type \c data_t
2115 
2116  This function is used, for example, in \ref o2scl::pinside.
2117 
2118  \note This function is not the same as a Givens rotation,
2119  which is typically referred to in BLAS routines as <tt>drot()</tt>.
2120  */
2121  template<class vec_t, class data_t>
2122  void vector_rotate(size_t n, vec_t &data, size_t k) {
2123 
2124  data_t *tmp=new data_t[n];
2125  for(size_t i=0;i<n;i++) {
2126  tmp[i]=data[(i+k)%n];
2127  }
2128  for(size_t i=0;i<n;i++) {
2129  data[i]=tmp[i];
2130  }
2131  delete[] tmp;
2132 
2133  return;
2134  }
2135 
2136  /** \brief Reverse the first \c n elements of a vector
2137 
2138  If \c n is zero, this function will silently do nothing.
2139  */
2140  template<class vec_t, class data_t>
2141  void vector_reverse(size_t n, vec_t &data) {
2142  data_t tmp;
2143 
2144  for(size_t i=0;i<n/2;i++) {
2145  tmp=data[n-1-i];
2146  data[n-1-i]=data[i];
2147  data[i]=tmp;
2148  }
2149  return;
2150  }
2151 
2152  /** \brief Reverse a vector
2153 
2154  If the <tt>size()</tt> method returns zero, this function will
2155  silently do nothing.
2156  */
2157  template<class vec_t, class data_t>
2158  void vector_reverse(vec_t &data) {
2159  data_t tmp;
2160  size_t n=data.size();
2161 
2162  for(size_t i=0;i<n/2;i++) {
2163  tmp=data[n-1-i];
2164  data[n-1-i]=data[i];
2165  data[i]=tmp;
2166  }
2167  return;
2168  }
2169 
2170  /** \brief Reverse the first n elements in a vector of double
2171  precision numbers
2172 
2173  If \c n is zero, this function will silently do nothing.
2174  */
2175  template<class vec_t>
2176  void vector_reverse_double(size_t n, vec_t &data) {
2177  double tmp;
2178 
2179  for(size_t i=0;i<n/2;i++) {
2180  tmp=data[n-1-i];
2181  data[n-1-i]=data[i];
2182  data[i]=tmp;
2183  }
2184  return;
2185  }
2186 
2187  /** \brief Reverse a vector of double precision numbers
2188 
2189  If the <tt>size()</tt> method returns zero, this function will
2190  silently do nothing.
2191  */
2192  template<class vec_t> void vector_reverse_double(vec_t &data) {
2193  double tmp;
2194  size_t n=data.size();
2195 
2196  for(size_t i=0;i<n/2;i++) {
2197  tmp=data[n-1-i];
2198  data[n-1-i]=data[i];
2199  data[i]=tmp;
2200  }
2201  return;
2202  }
2203 
2204  /** \brief Construct a row of a matrix
2205 
2206  This class template works with combinations of ublas
2207  <tt>matrix</tt> and <tt>matrix_row</tt> objects,
2208  <tt>arma::mat</tt> and <tt>arma::rowvec</tt>, and
2209  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2210 
2211  \note When calling this function with ublas objects, the
2212  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2213  otherwise some compilers will use argument dependent lookup and
2214  get (justifiably) confused with matrix_row in the ublas
2215  namespace.
2216 
2217  \note The template parameters must be explicitly specified
2218  when calling this template function.
2219  */
2220  template<class mat_t, class mat_row_t>
2221  mat_row_t matrix_row(mat_t &M, size_t row) {
2222  return mat_row_t(M,row);
2223  }
2224 
2225  /** \brief Generic object which represents a row of a matrix
2226 
2227  \note This class is experimental.
2228 
2229  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2230  to construct a row of a matrix object of type
2231  \code
2232  std::function<double &(size_t,size_t)>
2233  \endcode
2234  */
2235  template<class mat_t> class matrix_row_gen {
2236 
2237  protected:
2238 
2239  mat_t &m_;
2240 
2241  size_t row_;
2242 
2243  public:
2244 
2245  /// Create a row object from row \c row of matrix \c m
2246  matrix_row_gen(mat_t &m, size_t row) : m_(m), row_(row) {
2247  }
2248 
2249  /// Return a reference to the ith column of the selected row
2250  double &operator[](size_t i) {
2251  return m_(row_,i);
2252  }
2253 
2254  /// Return a const reference to the ith column of the selected row
2255  const double &operator[](size_t i) const {
2256  return m_(row_,i);
2257  }
2258  };
2259 
2260  /** \brief Construct a column of a matrix
2261 
2262  This class template works with combinations of ublas
2263  <tt>matrix</tt> and <tt>matrix_column</tt> objects,
2264  <tt>arma::mat</tt> and <tt>arma::colvec</tt>, and
2265  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2266 
2267  \note When calling this function with ublas objects, the
2268  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2269  otherwise some compilers will use argument dependent lookup and
2270  get (justifiably) confused with matrix_column in the ublas
2271  namespace.
2272 
2273  \note The template parameters must be explicitly specified
2274  when calling this template function.
2275  */
2276  template<class mat_t, class mat_column_t>
2277  mat_column_t matrix_column(mat_t &M, size_t column) {
2278  return mat_column_t(M,column);
2279  }
2280 
2281  /** \brief Generic object which represents a column of a matrix
2282 
2283  \note This class is experimental.
2284 
2285  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2286  to construct a row of a matrix object of type
2287  \code
2288  std::function<double &(size_t,size_t)>
2289  \endcode
2290  */
2291  template<class mat_t> class matrix_column_gen {
2292  protected:
2293  mat_t &m_;
2294  size_t column_;
2295  public:
2296  matrix_column_gen(mat_t &m, size_t column) : m_(m), column_(column) {
2297  }
2298  double &operator[](size_t i) {
2299  return m_(i,column_);
2300  }
2301  const double &operator[](size_t i) const {
2302  return m_(i,column_);
2303  }
2304  };
2305 
2306  /** \brief Output the first \c n elements of a vector to a stream,
2307  \c os
2308 
2309  No trailing space is output after the last element, and an
2310  endline is output only if \c endline is set to \c true. If the
2311  parameter \c n is zero, this function silently does nothing.
2312 
2313  This works with any class <tt>vec_t</tt> which has an operator[]
2314  which returns either the value of or a reference to the ith
2315  element and the element type has its own output operator which
2316  has been defined.
2317  */
2318  template<class vec_t>
2319  void vector_out(std::ostream &os, size_t n, const vec_t &v,
2320  bool endline=false) {
2321 
2322  // This next line is important since n-1 is not well-defined if n=0
2323  if (n==0) return;
2324 
2325  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2326  os << v[n-1];
2327  if (endline) os << std::endl;
2328  return;
2329  }
2330 
2331  /** \brief Output a vector to a stream
2332 
2333  No trailing space is output after the last element, and an
2334  endline is output only if \c endline is set to \c true. If the
2335  parameter \c n is zero, this function silently does nothing.
2336 
2337  This works with any class <tt>vec_t</tt> which has an operator[]
2338  which returns either the value of or a reference to the ith
2339  element and the element type has its own output operator which
2340  has been defined.
2341  */
2342  template<class vec_t>
2343  void vector_out(std::ostream &os, const vec_t &v, bool endline=false) {
2344 
2345  size_t n=v.size();
2346 
2347  // This next line is important since n-1 is not well-defined if n=0
2348  if (n==0) return;
2349 
2350  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2351  os << v[n-1];
2352  if (endline) os << std::endl;
2353  return;
2354  }
2355 
2356  /** \brief Fill a vector with a specified grid
2357  */
2358  template<class vec_t, class data_t>
2359  void vector_grid(uniform_grid<data_t> g, vec_t &v) {
2360  g.template vector<vec_t>(v);
2361  return;
2362  }
2363 
2364  /// Set a matrix to unity on the diagonal and zero otherwise
2365  template<class mat_t>
2366  void matrix_set_identity(size_t M, size_t N, mat_t &m) {
2367  for(size_t i=0;i<M;i++) {
2368  for(size_t j=0;j<N;j++) {
2369  if (i==j) m(i,j)=1.0;
2370  else m(i,j)=0.0;
2371  }
2372  }
2373  }
2374  //@}
2375 
2376  /// \name Vector range classes and functions
2377  //@{
2378  /** \brief Vector range function for pointers
2379 
2380  \note In this case, the return type is the same as the
2381  type of the first parameter.
2382  */
2383  template<class dat_t> dat_t *vector_range
2384  (dat_t *v, size_t start, size_t last) {
2385  return v+start;
2386  }
2387 
2388  /** \brief Vector range function for const pointers
2389 
2390  \note In this case, the return type is the same as the
2391  type of the first parameter.
2392  */
2393  template<class dat_t> const dat_t *const_vector_range
2394  (const dat_t *v, size_t start, size_t last) {
2395  return v+start;
2396  }
2397 
2398  /** \brief Vector range function template for ublas vectors
2399 
2400  The element with index \c start in the original vector
2401  will become the first argument in the new vector, and
2402  the new vector will have size <tt>last-start</tt> .
2403 
2404  \note In this case, the return type is not the same as the
2405  type of the first parameter.
2406  */
2407  template<class dat_t> boost::numeric::ublas::vector_range
2410  size_t start, size_t last) {
2413  (v,boost::numeric::ublas::range(start,last));
2414  }
2415 
2416  /** \brief Const vector range function template for ublas vectors
2417 
2418  The element with index \c start in the original vector
2419  will become the first argument in the new vector, and
2420  the new vector will have size <tt>last-start</tt> .
2421 
2422  \note In this case, the return type is not the same as the
2423  type of the first parameter.
2424  */
2425  template<class dat_t> const boost::numeric::ublas::vector_range
2428  size_t start, size_t last) {
2431  (v,boost::numeric::ublas::range(start,last));
2432  }
2433 
2434  /** \brief Const vector range function template for const ublas
2435  vectors
2436 
2437  The element with index \c start in the original vector
2438  will become the first argument in the new vector, and
2439  the new vector will have size <tt>last-start</tt> .
2440 
2441  \note In this case, the return type is not the same as the
2442  type of the first parameter.
2443  */
2444  template<class dat_t> const boost::numeric::ublas::vector_range
2447  size_t start, size_t last) {
2450  (v,boost::numeric::ublas::range(start,last));
2451  }
2452 
2453  /** \brief Vector range function template for ublas vector
2454  ranges of ublas vectors
2455 
2456  The element with index \c start in the original vector
2457  will become the first argument in the new vector, and
2458  the new vector will have size <tt>last-start</tt> .
2459 
2460  \note In this case, the return type is not the same as the
2461  type of the first parameter.
2462  */
2463  template<class dat_t>
2467  vector_range
2470  size_t start, size_t last) {
2474  (v,boost::numeric::ublas::range(start,last));
2475  }
2476 
2477  /** \brief Const vector range function template for ublas vector
2478  ranges of ublas vectors
2479 
2480  The element with index \c start in the original vector
2481  will become the first argument in the new vector, and
2482  the new vector will have size <tt>last-start</tt> .
2483 
2484  \note In this case, the return type is not the same as the
2485  type of the first parameter.
2486  */
2487  template<class dat_t>
2494  size_t start, size_t last) {
2498  (v,boost::numeric::ublas::range(start,last));
2499  }
2500 
2501  /** \brief Const vector range function template for const
2502  ublas vector ranges of ublas vectors
2503 
2504  The element with index \c start in the original vector
2505  will become the first argument in the new vector, and
2506  the new vector will have size <tt>last-start</tt> .
2507 
2508  \note In this case, the return type is not the same as the
2509  type of the first parameter.
2510  */
2511  template<class dat_t>
2518  size_t start, size_t last) {
2522  (v,boost::numeric::ublas::range(start,last));
2523  }
2524 
2525  /** \brief Const vector range function template for const
2526  ublas vector ranges of const ublas vectors
2527 
2528  The element with index \c start in the original vector
2529  will become the first argument in the new vector, and
2530  the new vector will have size <tt>last-start</tt> .
2531 
2532  \note In this case, the return type is not the same as the
2533  type of the first parameter.
2534  */
2535  template<class dat_t>
2542  size_t start, size_t last) {
2546  (v,boost::numeric::ublas::range(start,last));
2547  }
2548 
2549  // Forward definition for friendship
2550  template<class vec_t> class const_vector_range_gen;
2551 
2552  /** \brief Experimental vector range object
2553  */
2554  template<class vec_t> class vector_range_gen {
2555 
2556  protected:
2557 
2558  friend class const_vector_range_gen<vec_t>;
2559 
2560  /// A reference to the original vector
2561  vec_t &v_;
2562 
2563  /// The index offset
2564  size_t start_;
2565 
2566  /// The end() iterator
2567  size_t last_;
2568 
2569  public:
2570 
2571  /// Create an object starting with index \c start in vector \c v
2572  vector_range_gen(vec_t &v, size_t start, size_t last) : v_(v),
2573  start_(start), last_(last) {
2574 #if !O2SCL_NO_RANGE_CHECK
2575  if (last<start) {
2576  O2SCL_ERR2("End before beginning in vector_range_gen::",
2577  "vector_range_gen(vec_t,size_t,size_t)",
2579  }
2580 #endif
2581  }
2582 
2583  /// Create an object from a previously constructed range object
2584  vector_range_gen(const vector_range_gen &v2, size_t start,
2585  size_t last) : v_(v2.v_),
2586  start_(start+v2.start_), last_(last+v2.start_) {
2587 #if !O2SCL_NO_RANGE_CHECK
2588  if (last<start) {
2589  O2SCL_ERR2("End before beginning in vector_range_gen::",
2590  "vector_range_gen(vector_range_gen,size_t,size_t)",
2592  }
2593  if (last>v2.last_) {
2594  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2595  "vector_range_gen(vector_range_gen,size_t,size_t)",
2597  }
2598 #endif
2599  }
2600 
2601  /// Return the vector size
2602  size_t size() const {
2603  return last_-start_;
2604  }
2605 
2606  /// Return a reference ith element
2607  double &operator[](size_t i) {
2608 #if !O2SCL_NO_RANGE_CHECK
2609  if (i+start_>=last_) {
2610  O2SCL_ERR("Index out of range in vector_range_gen::operator[].",
2612  }
2613 #endif
2614  return v_[i+start_];
2615  }
2616 
2617  /// Return a const reference ith element
2618  const double &operator[](size_t i) const {
2619 #if !O2SCL_NO_RANGE_CHECK
2620  if (i+start_>=last_) {
2621  O2SCL_ERR2("Index out of range in ",
2622  "vector_range_gen::operator[] const.",o2scl::exc_einval);
2623  }
2624 #endif
2625  return v_[i+start_];
2626  }
2627  };
2628 
2629  /** \brief Experimental const vector range object
2630  */
2631  template<class vec_t> class const_vector_range_gen {
2632 
2633  protected:
2634 
2635  /// A reference to the original vector
2636  const vec_t &v_;
2637 
2638  /// The index offset
2639  size_t start_;
2640 
2641  /// The end() iterator
2642  size_t last_;
2643 
2644  public:
2645 
2646  /// Create an object starting with index \c start in vector \c v
2647  const_vector_range_gen(const vec_t &v, size_t start, size_t last) : v_(v),
2648  start_(start), last_(last) {
2649 #if !O2SCL_NO_RANGE_CHECK
2650  if (last<start) {
2651  O2SCL_ERR2("End before beginning in vector_range_gen::",
2652  "vector_range_gen(vec_t,size_t,size_t)",
2654  }
2655 #endif
2656  }
2657 
2658  /// Create an object from a previously constructed range object
2660  size_t last) : v_(v2.v_),
2661  start_(start+v2.start_), last_(last+v2.start_) {
2662 #if !O2SCL_NO_RANGE_CHECK
2663  if (last<start) {
2664  O2SCL_ERR2("End before beginning in vector_range_gen::",
2665  "vector_range_gen(vector_range_gen,size_t,size_t)",
2667  }
2668  if (last>v2.last_) {
2669  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2670  "vector_range_gen(vector_range_gen,size_t,size_t)",
2672  }
2673 #endif
2674  }
2675 
2676  /// Create an object from a previously constructed range object
2678  size_t last) : v_(v2.v_),
2679  start_(start+v2.start_), last_(last+v2.start_) {
2680 #if !O2SCL_NO_RANGE_CHECK
2681  if (last<start) {
2682  O2SCL_ERR2("End before beginning in vector_range_gen::",
2683  "vector_range_gen(vector_range_gen,size_t,size_t)",
2685  }
2686  if (last>v2.last_) {
2687  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2688  "vector_range_gen(vector_range_gen,size_t,size_t)",
2690  }
2691 #endif
2692  }
2693 
2694  /// Return the vector size
2695  size_t size() const {
2696  return last_-start_;
2697  }
2698 
2699  /// Return a const reference ith element
2700  const double &operator[](size_t i) const {
2701 #if !O2SCL_NO_RANGE_CHECK
2702  if (i+start_>=last_) {
2703  O2SCL_ERR2("Index out of range in ",
2704  "vector_range_gen::operator[] const.",o2scl::exc_einval);
2705  }
2706 #endif
2707  return v_[i+start_];
2708  }
2709  };
2710 
2711  /** \brief Create a \ref o2scl::vector_range_gen object
2712  from a <tt>std::vector</tt>
2713  */
2714  template<class data_t> vector_range_gen<std::vector<data_t> >
2715  vector_range(std::vector<data_t> &v, size_t start, size_t last) {
2716  return vector_range_gen<std::vector<data_t> >(v,start,last);
2717  }
2718 
2719  /** \brief Create a \ref o2scl::vector_range_gen object
2720  from a <tt>std::vector</tt>
2721  */
2722  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
2723  const_vector_range(const std::vector<data_t> &v, size_t start,
2724  size_t last) {
2725  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
2726  }
2727 
2728  /** \brief Create a \ref o2scl::vector_range_gen object
2729  from a <tt>std::vector</tt>
2730  */
2731  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
2732  const_vector_range(std::vector<data_t> &v, size_t start,
2733  size_t last) {
2734  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
2735  }
2736 
2737  /** \brief Recursively create a \ref o2scl::vector_range_gen object
2738  from a vector range
2739  */
2740  template<class vec_t> vector_range_gen<vec_t>
2741  vector_range(vector_range_gen<vec_t> &v, size_t start, size_t last) {
2742  return vector_range_gen<vec_t>(v,start,last);
2743  }
2744 
2745  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2746  object from a vector range
2747  */
2748  template<class vec_t> const const_vector_range_gen<vec_t>
2750  size_t start, size_t last) {
2751  return const_vector_range_gen<vec_t>(v,start,last);
2752  }
2753 
2754  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2755  object from a const vector range
2756  */
2757  template<class vec_t> const const_vector_range_gen<vec_t>
2759  size_t start, size_t last) {
2760  return const_vector_range_gen<vec_t>(v,start,last);
2761  }
2762 
2763  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2764  object from a const vector range
2765  */
2766  template<class vec_t> const const_vector_range_gen<vec_t>
2768  size_t start, size_t last) {
2769  return const_vector_range_gen<vec_t>(v,start,last);
2770  }
2771 
2772  /** \brief Vector range function template for <tt>std::vector</tt>
2773 
2774  The element with index \c start in the original vector
2775  will become the first argument in the new vector, and
2776  the new vector will have size <tt>last-start</tt> .
2777 
2778  \note In this case, the return type is the same as the
2779  type of the first parameter.
2780  \note Unlike the ublas and pointer cases, this forces
2781  a copy.
2782  */
2783  template<class dat_t> std::vector<dat_t>
2784  vector_range_copy(const std::vector<dat_t> &v, size_t start, size_t last) {
2785  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
2786  }
2787 
2788  /** \brief Const vector range function template for <tt>std::vector</tt>
2789 
2790  The element with index \c start in the original vector
2791  will become the first argument in the new vector, and
2792  the new vector will have size <tt>last-start</tt> .
2793 
2794  \note In this case, the return type is the same as the
2795  type of the first parameter.
2796  \note Unlike the ublas and pointer cases, this forces
2797  a copy.
2798  */
2799  template<class dat_t> const std::vector<dat_t>
2800  vector_range_copy(const std::vector<dat_t> &v, size_t start,
2801  size_t last) {
2802  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
2803  }
2804 
2805  //@}
2806 
2807 }
2808 
2809 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
2810 
2811 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
2812 #include <armadillo>
2813 namespace o2scl {
2814 
2815  /// \name Armadillo specializations
2816  //@{
2817  /// Armadillo version of \ref matrix_max()
2818  double matrix_max(const arma::mat &data);
2819 
2820  /// Armadillo version of \ref matrix_min()
2821  double matrix_min(const arma::mat &data);
2822 
2823  /// Armadillo version of \ref matrix_row()
2824  template<> arma::subview_row<double>
2825  matrix_row<arma::mat,arma::subview_row<double> >
2826  (arma::mat &M, size_t row);
2827 
2828  /// Armadillo version of \ref matrix_column()
2829  template<> arma::subview_col<double>
2830  matrix_column<arma::mat,arma::subview_col<double> >
2831  (arma::mat &M, size_t column);
2832  //@}
2833 
2834 }
2835 
2836 #endif
2837 
2838 #if defined (O2SCL_EIGEN) || defined (DOXYGEN)
2839 #include <eigen3/Eigen/Dense>
2840 
2841 namespace o2scl {
2842 
2843  /// \name Eigen specializations
2844  //@{
2845  /// Eigen version of \ref matrix_max()
2846  double matrix_max(const Eigen::MatrixXd &data);
2847 
2848  /// Eigen version of \ref matrix_min()
2849  double matrix_min(const Eigen::MatrixXd &data);
2850 
2851  /// Eigen version of \ref matrix_row()
2852  template<> Eigen::MatrixXd::RowXpr
2853  matrix_row<Eigen::MatrixXd,Eigen::MatrixXd::RowXpr>
2854  (Eigen::MatrixXd &M, size_t row);
2855 
2856  /// Eigen version of \ref matrix_column()
2857  template<> Eigen::MatrixXd::ColXpr
2858  matrix_column<Eigen::MatrixXd,Eigen::MatrixXd::ColXpr>
2859  (Eigen::MatrixXd &M, size_t column);
2860  //@}
2861 
2862 }
2863 
2864 #endif
2865 
2866 #else
2867 
2868 #include <o2scl/vector_special.h>
2869 
2870 // End of "#if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)"
2871 #endif
2872 
2873 #ifdef DOXYGEN
2874 /** \brief Placeholder documentation of some related Boost objects
2875  */
2876 namespace boost {
2877  /** \brief Documentation of Boost::numeric objects
2878  */
2879  namespace numeric {
2880  /** \brief Documentation of uBlas objects
2881  */
2882  namespace ublas {
2883  /** \brief The default vector type from uBlas
2884 
2885  The uBlas types aren't documented here, but the full documentation
2886  is available at
2887  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
2888 
2889  Internally in \o2, this is often typedef'd using
2890  \code
2891  typedef boost::numeric::ublas::vector<double> ubvector;
2892  typedef boost::numeric::ublas::vector<size_t> ubvector_size_t;
2893  typedef boost::numeric::ublas::vector<int> ubvector_int;
2894  \endcode
2895 
2896  This is documented in \ref vector.h .
2897  */
2898  template<class T, class A> class vector {
2899  };
2900  /** \brief The default matrix type from uBlas
2901 
2902  The uBlas types aren't documented here, but the full documentation
2903  is available at
2904  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
2905 
2906  Internally in \o2, this is often typedef'd using
2907  \code
2908  typedef boost::numeric::ublas::matrix<double> ubmatrix;
2909  typedef boost::numeric::ublas::matrix<size_t> ubmatrix_size_t;
2910  typedef boost::numeric::ublas::matrix<int> ubmatrix_int;
2911  \endcode
2912 
2913  This is documented in \ref vector.h .
2914  */
2915  template<class T, class F, class A> class matrix {
2916  };
2917  }
2918  }
2919 }
2920 // End of "#ifdef DOXYGEN"
2921 #endif
2922 
2923 // End of "#ifndef O2SCL_VECTOR_H"
2924 #endif
data_t vector_min_quad(size_t n, const vec_t &data)
Minimum of vector by quadratic fit.
Definition: vector.h:1344
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:1737
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:1189
Experimental vector range object.
Definition: vector.h:2554
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:1153
data_t matrix_sum(size_t m, size_t n, const mat_t &data)
Compute the sum of matrix elements.
Definition: vector.h:1702
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:1334
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:1369
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
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:626
void matrix_make_upper(mat_t &src)
Make a matrix upper triangular by setting the lower triangular entries to zero.
Definition: vector.h:349
Placeholder documentation of some related Boost objects.
Definition: vector.h:2876
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:1908
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:1451
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:2647
Generic object which represents a column of a matrix.
Definition: vector.h:2291
bool matrix_is_upper(mat_t &src)
Simple test that a matrix is upper triangular.
Definition: vector.h:320
bool matrix_is_lower(mat_t &src)
Simple test that a matrix is lower triangular.
Definition: vector.h:306
double & operator[](size_t i)
Return a reference ith element.
Definition: vector.h:2607
void vector_smallest_index(size_t n, 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:950
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:1863
A simple convenience wrapper for GSL vector objects.
Definition: vector.h:71
invalid argument supplied by user
Definition: err_hnd.h:59
Generic object which represents a row of a matrix.
Definition: vector.h:2235
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:2741
void vector_reverse(size_t n, vec_t &data)
Reverse the first n elements of a vector.
Definition: vector.h:2141
size_t last_
The end() iterator.
Definition: vector.h:2567
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:2122
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:493
size_t last_
The end() iterator.
Definition: vector.h:2642
double matrix_min_value_double(const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1554
The default matrix type from uBlas.
Definition: vector.h:2915
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:1208
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:2069
generic failure
Definition: err_hnd.h:61
void matrix_make_lower(mat_t &src)
Make a matrix lower triangular by setting the upper triangular entries to zero.
Definition: vector.h:335
double matrix_max_value_double(const mat_t &data)
Compute the maximum of a matrix.
Definition: vector.h:1427
data_t vector_sum(size_t n, vec_t &data)
Compute the sum of the first n elements of a vector.
Definition: vector.h:1968
void vector_swap(size_t N, vec_t &v1, vec2_t &v2)
Swap the first N elements of two vectors.
Definition: vector.h:423
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:1015
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:2784
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:2677
matrix_row_gen(mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2246
vec_t & v_
A reference to the original vector.
Definition: vector.h:2561
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:564
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts a vector (in increasing order)
Definition: vector.h:801
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:598
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
void matrix_copy(mat_t &src, mat2_t &dest)
Simple matrix copy.
Definition: vector.h:179
size_t size() const
Return the vector size.
Definition: vector.h:2695
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:1669
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:1078
const vec_t & v_
A reference to the original vector.
Definition: vector.h:2636
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2384
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
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:2659
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:1229
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:2366
data_t vector_max_quad(size_t n, const vec_t &data)
Maximum of vector by quadratic fit.
Definition: vector.h:1309
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:2584
mat_column_t matrix_column(mat_t &M, size_t column)
Construct a column of a matrix.
Definition: vector.h:2277
size_t start_
The index offset.
Definition: vector.h:2564
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
mat_row_t matrix_row(mat_t &M, size_t row)
Construct a row of a matrix.
Definition: vector.h:2221
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:2319
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:877
void vector_grid(uniform_grid< data_t > g, vec_t &v)
Fill a vector with a specified grid.
Definition: vector.h:2359
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:1279
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:2025
A simple convenience wrapper for GSL matrix objects.
Definition: vector.h:92
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:1114
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:1133
const double & operator[](size_t i) const
Return a const reference to the ith column of the selected row.
Definition: vector.h:2255
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:1252
double matrix_max(const arma::mat &data)
Armadillo version of matrix_max()
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:2618
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:2700
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:1996
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:2572
A class representing a uniform linear or logarithmic grid.
Definition: uniform_grid.h:38
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:1636
Experimental const vector range object.
Definition: vector.h:2550
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:1934
void matrix_transpose(mat_t &src, mat2_t &dest)
Simple transpose.
Definition: vector.h:225
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:737
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:2176
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:615
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:1579
double matrix_min(const arma::mat &data)
Armadillo version of matrix_min()
void matrix_swap(size_t M, size_t N, mat_t &v1, mat2_t &v2)
Swap of the first elements in two matrices.
Definition: vector.h:545
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:643
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:1382
size_t start_
The index offset.
Definition: vector.h:2639
std::string itos(int x)
Convert an integer to a string.
std::string szttos(size_t x)
Convert a size_t to a string.
const dat_t * const_vector_range(const dat_t *v, size_t start, size_t last)
Vector range function for const pointers.
Definition: vector.h:2394
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:1784
double & operator[](size_t i)
Return a reference to the ith column of the selected row.
Definition: vector.h:2250
data_t matrix_min_value(size_t m, size_t n, const mat_t &data)
Compute the minimum of a matrix.
Definition: vector.h:1508
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Definition: vector.h:852
The default vector type from uBlas.
Definition: vector.h:2898
size_t size() const
Return the vector size.
Definition: vector.h:2602
void vector_sort(size_t n, vec_t &data)
Sort a vector (in increasing order)
Definition: vector.h:708
void sort_downheap(vec_t &data, size_t n, size_t k)
Provide a downheap() function for vector_sort()
Definition: vector.h:653

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