vector.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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 Ordering and finite tests
1961  //@{
1962  /** \brief Test if the first \c n elements of a vector are
1963  monotonic and increasing or decreasing
1964 
1965  If \c n is zero or one, this function will return 0 without
1966  calling the error handler. If all the vector's elements are equal,
1967  this function will return 3. Otherwise, if the vector is not
1968  monotonic, then this function will return 0. Finally, if the
1969  vector is nondecreasing (increasing or equal intervals), this
1970  function will return 1, and if the vector is nonincreasing
1971  (decreasing or equal intervals), this function will return 2.
1972  This function assumes that simple comparison operators have been
1973  defined for the type of each vector element.
1974  */
1975  template<class vec_t>
1976  int vector_is_monotonic(size_t n, vec_t &data) {
1977 
1978  if (n<1) return 0;
1979  if (n<2) {
1980  if (data[0]==data[1]) {
1981  return 3;
1982  } else if (data[0]<data[1]) {
1983  return 1;
1984  } else {
1985  return 2;
1986  }
1987  }
1988 
1989  // Find first non-flat interval
1990  size_t start=0;
1991  bool done=false;
1992  for(size_t i=0;i<n-1 && done==false;i++) {
1993  if (data[i]!=data[i+1]) {
1994  done=true;
1995  } else {
1996  start++;
1997  }
1998  }
1999 
2000  // If all elements in the vector are equal, or if the only
2001  // distinct element is at the end, then return true.
2002  if (done==false) {
2003  return 3;
2004  }
2005 
2006  if (start==n-2) {
2007  if (data[start]<data[start+1]) return 1;
2008  else return 2;
2009  }
2010 
2011  // Determine if the vector is increasing (inc=true) or decreasing
2012  // (inc=false)
2013  bool inc=true;
2014  if (data[start]>data[start+1]) inc=false;
2015 
2016  if (inc) {
2017  for(size_t i=start+1;i<n-1;i++) {
2018  // If there is one decreasing interval, return false
2019  if (data[i]>data[i+1]) return 0;
2020  }
2021  return 1;
2022  }
2023 
2024  // If there is one increasing interval, return false
2025  for(size_t i=start+1;i<n-1;i++) {
2026  if (data[i]<data[i+1]) return 0;
2027  }
2028  return 2;
2029  }
2030 
2031  /** \brief Test if the first \c n elements of a vector are
2032  monotonic and increasing or decreasing
2033 
2034  If \c n is zero or one, this function will return 0 without
2035  calling the error handler. If all the vector's elements are equal,
2036  this function will return 3. Otherwise, if the vector is not
2037  monotonic, then this function will return 0. Finally, if the
2038  vector is nondecreasing (increasing or equal intervals), this
2039  function will return 1, and if the vector is nonincreasing
2040  (decreasing or equal intervals), this function will return 2.
2041  This function assumes that simple comparison operators have been
2042  defined for the type of each vector element.
2043  */
2044  template<class vec_t> int vector_is_monotonic(vec_t &data) {
2045  return vector_is_monotonic(data.size(),data);
2046  }
2047 
2048  /** \brief Test if the first \c n elements of a vector are
2049  strictly monotonic and determine if they are increasing or decreasing
2050 
2051  If \c n is zero this function will return 0 without calling the
2052  error handler. Also, if the vector is not monotonic, this
2053  function will return 0. If the vector is strictly
2054  monotonic, then this function will return 1 if it is
2055  increasing and 2 if it is decreasing.
2056  */
2057  template<class vec_t>
2058  int vector_is_strictly_monotonic(size_t n, vec_t &data) {
2059 
2060  if (n<1) return 0;
2061 
2062  // Determine if the vector is increasing (inc=true) or decreasing
2063  // (inc=false)
2064  bool inc=true;
2065  if (data[0]==data[1]) {
2066  return 0;
2067  } else if (data[0]>data[1]) {
2068  inc=false;
2069  }
2070 
2071  if (inc) {
2072  for(size_t i=1;i<n-1;i++) {
2073  // If there is one nonincreasing interval, return 0
2074  if (data[i]>=data[i+1]) return 0;
2075  }
2076  return 1;
2077  }
2078 
2079  // If there is one increasing interval, return 0
2080  for(size_t i=1;i<n-1;i++) {
2081  if (data[i]<=data[i+1]) return 0;
2082  }
2083  return 2;
2084  }
2085 
2086  /** \brief Test if the first \c n elements of a vector are
2087  strictly monotonic and determine if they are increasing or decreasing
2088 
2089  If \c n is zero this function will return 0 without calling the
2090  error handler. Also, if the vector is not monotonic, this
2091  function will return 0. If the vector is strictly
2092  monotonic, then this function will return 1 if it is
2093  increasing and 2 if it is decreasing.
2094  */
2095  template<class vec_t>
2096  int vector_is_strictly_monotonic(vec_t &data) {
2097  return vector_is_strictly_monotonic(data.size(),data);
2098  }
2099 
2100  /** \brief Test if the first \c n elements of a vector are finite
2101 
2102  If \c n is zero, this will return true without throwing
2103  an exception.
2104 
2105  The corresponding tests for matrix functions are
2106  in clbas_base.h .
2107  */
2108  template<class vec_t>
2109  bool vector_is_finite(size_t n, vec_t &data) {
2110  for(size_t i=0;i<n;i++) {
2111  if (!std::isfinite(data[i])) return false;
2112  }
2113  return true;
2114  }
2115 
2116  /** \brief Test if a vector is finite
2117 
2118  If \c n is zero, this will return true without throwing
2119  an exception.
2120 
2121  The corresponding tests for matrix functions are
2122  in clbas_base.h .
2123  */
2124  template<class vec_t> bool vector_is_finite(vec_t &data) {
2125  return vector_is_finite(data.size(),data);
2126  }
2127  //@}
2128 
2129  /// \name Miscellaneous mathematical functions
2130  //@{
2131  /** \brief Compute the sum of the first \c n elements of a vector
2132 
2133  If \c n is zero, this will return 0 without throwing
2134  an exception.
2135  */
2136  template<class vec_t, class data_t>
2137  data_t vector_sum(size_t n, vec_t &data) {
2138 
2139  data_t sum=0.0;
2140  for(size_t i=0;i<n;i++) {
2141  sum+=data[i];
2142  }
2143  return sum;
2144  }
2145 
2146  /** \brief Compute the sum of all the elements of a vector
2147 
2148  If the vector has zero size, this will return 0 without
2149  calling the error handler.
2150  */
2151  template<class vec_t, class data_t> data_t vector_sum(vec_t &data) {
2152  data_t sum=0.0;
2153  for(size_t i=0;i<data.size();i++) {
2154  sum+=data[i];
2155  }
2156  return sum;
2157  }
2158 
2159  /** \brief Compute the sum of the first \c n elements of a vector
2160  of double-precision numbers
2161 
2162  If \c n is zero, this will return 0 without throwing
2163  an exception.
2164  */
2165  template<class vec_t>double vector_sum_double(size_t n, vec_t &data) {
2166  double sum=0.0;
2167  for(size_t i=0;i<n;i++) {
2168  sum+=data[i];
2169  }
2170  return sum;
2171  }
2172 
2173  /** \brief Compute the sum of all the elements of a vector
2174  of double-precision numbers
2175 
2176  If the vector has zero size, this will return 0 without
2177  calling the error handler.
2178  */
2179  template<class vec_t> double vector_sum_double(vec_t &data) {
2180  double sum=0.0;
2181  for(size_t i=0;i<data.size();i++) {
2182  sum+=data[i];
2183  }
2184  return sum;
2185  }
2186 
2187  /** \brief Compute the norm of the first \c n entries of a
2188  vector of floating-point (single or double precision) numbers
2189 
2190  This function is a more generic version of
2191  \ref o2scl_cblas::dnrm2 .
2192  */
2193  template<class vec_t, class data_t>
2194  data_t vector_norm(size_t n, const vec_t &x) {
2195 
2196  data_t scale = 0.0;
2197  data_t ssq = 1.0;
2198 
2199  if (n <= 0) {
2200  return 0.0;
2201  } else if (n == 1) {
2202  return fabs(x[0]);
2203  }
2204 
2205  for (size_t i = 0; i < n; i++) {
2206  const data_t xx = x[i];
2207 
2208  if (xx != 0.0) {
2209  const data_t ax = fabs(xx);
2210 
2211  if (scale < ax) {
2212  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2213  scale = ax;
2214  } else {
2215  ssq += (ax / scale) * (ax / scale);
2216  }
2217  }
2218 
2219  }
2220 
2221  return scale * sqrt(ssq);
2222  }
2223 
2224  /** \brief Compute the norm of a vector of floating-point
2225  (single or double precision) numbers
2226  */
2227  template<class vec_t, class data_t> data_t vector_norm(const vec_t &x) {
2228  return vector_norm<vec_t,data_t>(x.size(),x);
2229  }
2230 
2231  /** \brief Compute the norm of the first \c n entries of a
2232  vector of double precision numbers
2233 
2234  This function is a more generic version of
2235  \ref o2scl_cblas::dnrm2 .
2236  */
2237  template<class vec_t>
2238  double vector_norm_double(size_t n, const vec_t &x) {
2239 
2240  double scale = 0.0;
2241  double ssq = 1.0;
2242 
2243  if (n <= 0) {
2244  return 0.0;
2245  } else if (n == 1) {
2246  return fabs(x[0]);
2247  }
2248 
2249  for (size_t i = 0; i < n; i++) {
2250  const double xx = x[i];
2251 
2252  if (xx != 0.0) {
2253  const double ax = fabs(xx);
2254 
2255  if (scale < ax) {
2256  ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
2257  scale = ax;
2258  } else {
2259  ssq += (ax / scale) * (ax / scale);
2260  }
2261  }
2262 
2263  }
2264 
2265  return scale * sqrt(ssq);
2266  }
2267 
2268  /** \brief Compute the norm of a vector of double precision numbers
2269  */
2270  template<class vec_t> double vector_norm_double(const vec_t &x) {
2271  return vector_norm_double<vec_t>(x.size(),x);
2272  }
2273  //@}
2274 
2275  /// \name Other vector and matrix functions
2276  //@{
2277  /** \brief Set the first N entries in a vector to a particular value
2278  */
2279  template<class vec_t, class data_t>
2280  void vector_set_all(size_t N, vec_t &src, data_t &val) {
2281  for(size_t i=0;i<N;i++) {
2282  src[i]=val;
2283  }
2284  return;
2285  }
2286 
2287  /** \brief Set all entries in a vector to a particular value
2288  */
2289  template<class vec_t, class data_t>
2290  void vector_set_all(vec_t &src, data_t &val) {
2291  o2scl::vector_set_all(src.size(),src,val);
2292  return;
2293  }
2294 
2295  /** \brief From a given vector, create a new vector by removing a
2296  specified element
2297  */
2298  template<class vec_t, class vec2_t>
2299  void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest) {
2300  if (src.size()==0) {
2301  O2SCL_ERR("Empty source vector.",o2scl::exc_einval);
2302  }
2303  if (iout>=src.size()) {
2304  O2SCL_ERR("Requested element beyond end.",o2scl::exc_einval);
2305  }
2306  dest.resize(src.size()-1);
2307  size_t j=0;
2308  for(size_t i=0;i<src.size();i++) {
2309  if (i!=iout) {
2310  dest[j]=src[i];
2311  j++;
2312  }
2313  }
2314  return;
2315  }
2316 
2317  /** \brief "Rotate" a vector so that the kth element is now the beginning
2318 
2319  This is a generic template function which will work for
2320  any types \c data_t and \c vec_t for which
2321  - \c data_t has an <tt>operator=</tt>
2322  - <tt>vec_t::operator[]</tt> returns a reference
2323  to an object of type \c data_t
2324 
2325  This function is used, for example, in \ref o2scl::pinside.
2326 
2327  \note This function is not the same as a Givens rotation,
2328  which is typically referred to in BLAS routines as <tt>drot()</tt>.
2329  */
2330  template<class vec_t, class data_t>
2331  void vector_rotate(size_t n, vec_t &data, size_t k) {
2332 
2333  data_t *tmp=new data_t[n];
2334  for(size_t i=0;i<n;i++) {
2335  tmp[i]=data[(i+k)%n];
2336  }
2337  for(size_t i=0;i<n;i++) {
2338  data[i]=tmp[i];
2339  }
2340  delete[] tmp;
2341 
2342  return;
2343  }
2344 
2345  /** \brief Reverse the first \c n elements of a vector
2346 
2347  If \c n is zero, this function will silently do nothing.
2348  */
2349  template<class vec_t, class data_t>
2350  void vector_reverse(size_t n, vec_t &data) {
2351  data_t tmp;
2352 
2353  for(size_t i=0;i<n/2;i++) {
2354  tmp=data[n-1-i];
2355  data[n-1-i]=data[i];
2356  data[i]=tmp;
2357  }
2358  return;
2359  }
2360 
2361  /** \brief Reverse a vector
2362 
2363  If the <tt>size()</tt> method returns zero, this function will
2364  silently do nothing.
2365  */
2366  template<class vec_t, class data_t>
2367  void vector_reverse(vec_t &data) {
2368  data_t tmp;
2369  size_t n=data.size();
2370 
2371  for(size_t i=0;i<n/2;i++) {
2372  tmp=data[n-1-i];
2373  data[n-1-i]=data[i];
2374  data[i]=tmp;
2375  }
2376  return;
2377  }
2378 
2379  /** \brief Reverse the first n elements in a vector of double
2380  precision numbers
2381 
2382  If \c n is zero, this function will silently do nothing.
2383  */
2384  template<class vec_t>
2385  void vector_reverse_double(size_t n, vec_t &data) {
2386  double tmp;
2387 
2388  for(size_t i=0;i<n/2;i++) {
2389  tmp=data[n-1-i];
2390  data[n-1-i]=data[i];
2391  data[i]=tmp;
2392  }
2393  return;
2394  }
2395 
2396  /** \brief Reverse a vector of double precision numbers
2397 
2398  If the <tt>size()</tt> method returns zero, this function will
2399  silently do nothing.
2400  */
2401  template<class vec_t> void vector_reverse_double(vec_t &data) {
2402  double tmp;
2403  size_t n=data.size();
2404 
2405  for(size_t i=0;i<n/2;i++) {
2406  tmp=data[n-1-i];
2407  data[n-1-i]=data[i];
2408  data[i]=tmp;
2409  }
2410  return;
2411  }
2412 
2413  /** \brief Construct a row of a matrix
2414 
2415  This class template works with combinations of ublas
2416  <tt>matrix</tt> and <tt>matrix_row</tt> objects,
2417  <tt>arma::mat</tt> and <tt>arma::rowvec</tt>, and
2418  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2419 
2420  \note When calling this function with ublas objects, the
2421  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2422  otherwise some compilers will use argument dependent lookup and
2423  get (justifiably) confused with matrix_row in the ublas
2424  namespace.
2425 
2426  \note The template parameters must be explicitly specified
2427  when calling this template function.
2428  */
2429  template<class mat_t, class mat_row_t>
2430  mat_row_t matrix_row(mat_t &M, size_t row) {
2431  return mat_row_t(M,row);
2432  }
2433 
2434  /** \brief Generic object which represents a row of a matrix
2435 
2436  \note This class is experimental.
2437 
2438  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2439  to construct a row of a matrix object of type
2440  \code
2441  std::function<double &(size_t,size_t)>
2442  \endcode
2443  */
2444  template<class mat_t> class matrix_row_gen {
2445 
2446  protected:
2447 
2448  mat_t &m_;
2449 
2450  size_t row_;
2451 
2452  public:
2453 
2454  /// Create a row object from row \c row of matrix \c m
2455  matrix_row_gen(mat_t &m, size_t row) : m_(m), row_(row) {
2456  }
2457 
2458  /// Return a reference to the ith column of the selected row
2459  double &operator[](size_t i) {
2460  return m_(row_,i);
2461  }
2462 
2463  /// Return a const reference to the ith column of the selected row
2464  const double &operator[](size_t i) const {
2465  return m_(row_,i);
2466  }
2467  };
2468 
2469  /** \brief Construct a column of a matrix
2470 
2471  This class template works with combinations of ublas
2472  <tt>matrix</tt> and <tt>matrix_column</tt> objects,
2473  <tt>arma::mat</tt> and <tt>arma::colvec</tt>, and
2474  <tt>Eigen::MatrixXd</tt> and <tt>Eigen::VectorXd</tt>.
2475 
2476  \note When calling this function with ublas objects, the
2477  namespace prefix <tt>"o2scl::"</tt> must often be specified,
2478  otherwise some compilers will use argument dependent lookup and
2479  get (justifiably) confused with matrix_column in the ublas
2480  namespace.
2481 
2482  \note The template parameters must be explicitly specified
2483  when calling this template function.
2484  */
2485  template<class mat_t, class mat_column_t>
2486  mat_column_t matrix_column(mat_t &M, size_t column) {
2487  return mat_column_t(M,column);
2488  }
2489 
2490  /** \brief Generic object which represents a column of a matrix
2491 
2492  \note This class is experimental.
2493 
2494  This class is used in <tt>o2scl::eos_sn_base::slice</tt>
2495  to construct a row of a matrix object of type
2496  \code
2497  std::function<double &(size_t,size_t)>
2498  \endcode
2499  */
2500  template<class mat_t> class matrix_column_gen {
2501  protected:
2502  mat_t &m_;
2503  size_t column_;
2504  public:
2505  matrix_column_gen(mat_t &m, size_t column) : m_(m), column_(column) {
2506  }
2507  double &operator[](size_t i) {
2508  return m_(i,column_);
2509  }
2510  const double &operator[](size_t i) const {
2511  return m_(i,column_);
2512  }
2513  };
2514 
2515  /** \brief Output the first \c n elements of a vector to a stream,
2516  \c os
2517 
2518  No trailing space is output after the last element, and an
2519  endline is output only if \c endline is set to \c true. If the
2520  parameter \c n is zero, this function silently does nothing.
2521 
2522  This works with any class <tt>vec_t</tt> which has an operator[]
2523  which returns either the value of or a reference to the ith
2524  element and the element type has its own output operator which
2525  has been defined.
2526  */
2527  template<class vec_t>
2528  void vector_out(std::ostream &os, size_t n, const vec_t &v,
2529  bool endline=false) {
2530 
2531  // This next line is important since n-1 is not well-defined if n=0
2532  if (n==0) return;
2533 
2534  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2535  os << v[n-1];
2536  if (endline) os << std::endl;
2537  return;
2538  }
2539 
2540  /** \brief Output a vector to a stream
2541 
2542  No trailing space is output after the last element, and an
2543  endline is output only if \c endline is set to \c true. If the
2544  parameter \c n is zero, this function silently does nothing.
2545 
2546  This works with any class <tt>vec_t</tt> which has an operator[]
2547  which returns either the value of or a reference to the ith
2548  element and the element type has its own output operator which
2549  has been defined.
2550  */
2551  template<class vec_t>
2552  void vector_out(std::ostream &os, const vec_t &v, bool endline=false) {
2553 
2554  size_t n=v.size();
2555 
2556  // This next line is important since n-1 is not well-defined if n=0
2557  if (n==0) return;
2558 
2559  for(size_t i=0;i<n-1;i++) os << v[i] << " ";
2560  os << v[n-1];
2561  if (endline) os << std::endl;
2562  return;
2563  }
2564 
2565  /** \brief Fill a vector with a specified grid
2566  */
2567  template<class vec_t, class data_t>
2568  void vector_grid(uniform_grid<data_t> g, vec_t &v) {
2569  g.template vector<vec_t>(v);
2570  return;
2571  }
2572 
2573  /// Set a matrix to unity on the diagonal and zero otherwise
2574  template<class mat_t>
2575  void matrix_set_identity(size_t M, size_t N, mat_t &m) {
2576  for(size_t i=0;i<M;i++) {
2577  for(size_t j=0;j<N;j++) {
2578  if (i==j) m(i,j)=1.0;
2579  else m(i,j)=0.0;
2580  }
2581  }
2582  }
2583  //@}
2584 
2585  /// \name Vector range classes and functions
2586  //@{
2587  /** \brief Vector range function for pointers
2588 
2589  \note In this case, the return type is the same as the
2590  type of the first parameter.
2591  */
2592  template<class dat_t> dat_t *vector_range
2593  (dat_t *v, size_t start, size_t last) {
2594  return v+start;
2595  }
2596 
2597  /** \brief Vector range function for const pointers
2598 
2599  \note In this case, the return type is the same as the
2600  type of the first parameter.
2601  */
2602  template<class dat_t> const dat_t *const_vector_range
2603  (const dat_t *v, size_t start, size_t last) {
2604  return v+start;
2605  }
2606 
2607  /** \brief Vector range function template for ublas vectors
2608 
2609  The element with index \c start in the original vector
2610  will become the first argument in the new vector, and
2611  the new vector will have size <tt>last-start</tt> .
2612 
2613  \note In this case, the return type is not the same as the
2614  type of the first parameter.
2615  */
2616  template<class dat_t> boost::numeric::ublas::vector_range
2619  size_t start, size_t last) {
2622  (v,boost::numeric::ublas::range(start,last));
2623  }
2624 
2625  /** \brief Const vector range function template for ublas vectors
2626 
2627  The element with index \c start in the original vector
2628  will become the first argument in the new vector, and
2629  the new vector will have size <tt>last-start</tt> .
2630 
2631  \note In this case, the return type is not the same as the
2632  type of the first parameter.
2633  */
2634  template<class dat_t> const boost::numeric::ublas::vector_range
2637  size_t start, size_t last) {
2640  (v,boost::numeric::ublas::range(start,last));
2641  }
2642 
2643  /** \brief Const vector range function template for const ublas
2644  vectors
2645 
2646  The element with index \c start in the original vector
2647  will become the first argument in the new vector, and
2648  the new vector will have size <tt>last-start</tt> .
2649 
2650  \note In this case, the return type is not the same as the
2651  type of the first parameter.
2652  */
2653  template<class dat_t> const boost::numeric::ublas::vector_range
2656  size_t start, size_t last) {
2659  (v,boost::numeric::ublas::range(start,last));
2660  }
2661 
2662  /** \brief Vector range function template for ublas vector
2663  ranges of ublas vectors
2664 
2665  The element with index \c start in the original vector
2666  will become the first argument in the new vector, and
2667  the new vector will have size <tt>last-start</tt> .
2668 
2669  \note In this case, the return type is not the same as the
2670  type of the first parameter.
2671  */
2672  template<class dat_t>
2676  vector_range
2679  size_t start, size_t last) {
2683  (v,boost::numeric::ublas::range(start,last));
2684  }
2685 
2686  /** \brief Const vector range function template for ublas vector
2687  ranges of ublas vectors
2688 
2689  The element with index \c start in the original vector
2690  will become the first argument in the new vector, and
2691  the new vector will have size <tt>last-start</tt> .
2692 
2693  \note In this case, the return type is not the same as the
2694  type of the first parameter.
2695  */
2696  template<class dat_t>
2703  size_t start, size_t last) {
2707  (v,boost::numeric::ublas::range(start,last));
2708  }
2709 
2710  /** \brief Const vector range function template for const
2711  ublas vector ranges of ublas vectors
2712 
2713  The element with index \c start in the original vector
2714  will become the first argument in the new vector, and
2715  the new vector will have size <tt>last-start</tt> .
2716 
2717  \note In this case, the return type is not the same as the
2718  type of the first parameter.
2719  */
2720  template<class dat_t>
2727  size_t start, size_t last) {
2731  (v,boost::numeric::ublas::range(start,last));
2732  }
2733 
2734  /** \brief Const vector range function template for const
2735  ublas vector ranges of const ublas vectors
2736 
2737  The element with index \c start in the original vector
2738  will become the first argument in the new vector, and
2739  the new vector will have size <tt>last-start</tt> .
2740 
2741  \note In this case, the return type is not the same as the
2742  type of the first parameter.
2743  */
2744  template<class dat_t>
2751  size_t start, size_t last) {
2755  (v,boost::numeric::ublas::range(start,last));
2756  }
2757 
2758  // Forward definition for friendship
2759  template<class vec_t> class const_vector_range_gen;
2760 
2761  /** \brief Experimental vector range object
2762  */
2763  template<class vec_t> class vector_range_gen {
2764 
2765  protected:
2766 
2767  friend class const_vector_range_gen<vec_t>;
2768 
2769  /// A reference to the original vector
2770  vec_t &v_;
2771 
2772  /// The index offset
2773  size_t start_;
2774 
2775  /// The end() iterator
2776  size_t last_;
2777 
2778  public:
2779 
2780  /// Create an object starting with index \c start in vector \c v
2781  vector_range_gen(vec_t &v, size_t start, size_t last) : v_(v),
2782  start_(start), last_(last) {
2783 #if !O2SCL_NO_RANGE_CHECK
2784  if (last<start) {
2785  O2SCL_ERR2("End before beginning in vector_range_gen::",
2786  "vector_range_gen(vec_t,size_t,size_t)",
2788  }
2789 #endif
2790  }
2791 
2792  /// Create an object from a previously constructed range object
2793  vector_range_gen(const vector_range_gen &v2, size_t start,
2794  size_t last) : v_(v2.v_),
2795  start_(start+v2.start_), last_(last+v2.start_) {
2796 #if !O2SCL_NO_RANGE_CHECK
2797  if (last<start) {
2798  O2SCL_ERR2("End before beginning in vector_range_gen::",
2799  "vector_range_gen(vector_range_gen,size_t,size_t)",
2801  }
2802  if (last>v2.last_) {
2803  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2804  "vector_range_gen(vector_range_gen,size_t,size_t)",
2806  }
2807 #endif
2808  }
2809 
2810  /// Return the vector size
2811  size_t size() const {
2812  return last_-start_;
2813  }
2814 
2815  /// Return a reference ith element
2816  double &operator[](size_t i) {
2817 #if !O2SCL_NO_RANGE_CHECK
2818  if (i+start_>=last_) {
2819  O2SCL_ERR("Index out of range in vector_range_gen::operator[].",
2821  }
2822 #endif
2823  return v_[i+start_];
2824  }
2825 
2826  /// Return a const reference ith element
2827  const double &operator[](size_t i) const {
2828 #if !O2SCL_NO_RANGE_CHECK
2829  if (i+start_>=last_) {
2830  O2SCL_ERR2("Index out of range in ",
2831  "vector_range_gen::operator[] const.",o2scl::exc_einval);
2832  }
2833 #endif
2834  return v_[i+start_];
2835  }
2836  };
2837 
2838  /** \brief Experimental const vector range object
2839  */
2840  template<class vec_t> class const_vector_range_gen {
2841 
2842  protected:
2843 
2844  /// A reference to the original vector
2845  const vec_t &v_;
2846 
2847  /// The index offset
2848  size_t start_;
2849 
2850  /// The end() iterator
2851  size_t last_;
2852 
2853  public:
2854 
2855  /// Create an object starting with index \c start in vector \c v
2856  const_vector_range_gen(const vec_t &v, size_t start, size_t last) : v_(v),
2857  start_(start), last_(last) {
2858 #if !O2SCL_NO_RANGE_CHECK
2859  if (last<start) {
2860  O2SCL_ERR2("End before beginning in vector_range_gen::",
2861  "vector_range_gen(vec_t,size_t,size_t)",
2863  }
2864 #endif
2865  }
2866 
2867  /// Create an object from a previously constructed range object
2869  size_t last) : v_(v2.v_),
2870  start_(start+v2.start_), last_(last+v2.start_) {
2871 #if !O2SCL_NO_RANGE_CHECK
2872  if (last<start) {
2873  O2SCL_ERR2("End before beginning in vector_range_gen::",
2874  "vector_range_gen(vector_range_gen,size_t,size_t)",
2876  }
2877  if (last>v2.last_) {
2878  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2879  "vector_range_gen(vector_range_gen,size_t,size_t)",
2881  }
2882 #endif
2883  }
2884 
2885  /// Create an object from a previously constructed range object
2887  size_t last) : v_(v2.v_),
2888  start_(start+v2.start_), last_(last+v2.start_) {
2889 #if !O2SCL_NO_RANGE_CHECK
2890  if (last<start) {
2891  O2SCL_ERR2("End before beginning in vector_range_gen::",
2892  "vector_range_gen(vector_range_gen,size_t,size_t)",
2894  }
2895  if (last>v2.last_) {
2896  O2SCL_ERR2("End beyond end of previous vector in vector_range_gen::",
2897  "vector_range_gen(vector_range_gen,size_t,size_t)",
2899  }
2900 #endif
2901  }
2902 
2903  /// Return the vector size
2904  size_t size() const {
2905  return last_-start_;
2906  }
2907 
2908  /// Return a const reference ith element
2909  const double &operator[](size_t i) const {
2910 #if !O2SCL_NO_RANGE_CHECK
2911  if (i+start_>=last_) {
2912  O2SCL_ERR2("Index out of range in ",
2913  "vector_range_gen::operator[] const.",o2scl::exc_einval);
2914  }
2915 #endif
2916  return v_[i+start_];
2917  }
2918  };
2919 
2920  /** \brief Create a \ref o2scl::vector_range_gen object
2921  from a <tt>std::vector</tt>
2922  */
2923  template<class data_t> vector_range_gen<std::vector<data_t> >
2924  vector_range(std::vector<data_t> &v, size_t start, size_t last) {
2925  return vector_range_gen<std::vector<data_t> >(v,start,last);
2926  }
2927 
2928  /** \brief Create a \ref o2scl::vector_range_gen object
2929  from a <tt>std::vector</tt>
2930  */
2931  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
2932  const_vector_range(const std::vector<data_t> &v, size_t start,
2933  size_t last) {
2934  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
2935  }
2936 
2937  /** \brief Create a \ref o2scl::vector_range_gen object
2938  from a <tt>std::vector</tt>
2939  */
2940  template<class data_t> const const_vector_range_gen<std::vector<data_t> >
2941  const_vector_range(std::vector<data_t> &v, size_t start,
2942  size_t last) {
2943  return const_vector_range_gen<std::vector<data_t> >(v,start,last);
2944  }
2945 
2946  /** \brief Recursively create a \ref o2scl::vector_range_gen object
2947  from a vector range
2948  */
2949  template<class vec_t> vector_range_gen<vec_t>
2950  vector_range(vector_range_gen<vec_t> &v, size_t start, size_t last) {
2951  return vector_range_gen<vec_t>(v,start,last);
2952  }
2953 
2954  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2955  object from a vector range
2956  */
2957  template<class vec_t> const const_vector_range_gen<vec_t>
2959  size_t start, size_t last) {
2960  return const_vector_range_gen<vec_t>(v,start,last);
2961  }
2962 
2963  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2964  object from a const vector range
2965  */
2966  template<class vec_t> const const_vector_range_gen<vec_t>
2968  size_t start, size_t last) {
2969  return const_vector_range_gen<vec_t>(v,start,last);
2970  }
2971 
2972  /** \brief Recursively create a const \ref o2scl::vector_range_gen
2973  object from a const vector range
2974  */
2975  template<class vec_t> const const_vector_range_gen<vec_t>
2977  size_t start, size_t last) {
2978  return const_vector_range_gen<vec_t>(v,start,last);
2979  }
2980 
2981  /** \brief Vector range function template for <tt>std::vector</tt>
2982 
2983  The element with index \c start in the original vector
2984  will become the first argument in the new vector, and
2985  the new vector will have size <tt>last-start</tt> .
2986 
2987  \note In this case, the return type is the same as the
2988  type of the first parameter.
2989  \note Unlike the ublas and pointer cases, this forces
2990  a copy.
2991  */
2992  template<class dat_t> std::vector<dat_t>
2993  vector_range_copy(const std::vector<dat_t> &v, size_t start, size_t last) {
2994  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
2995  }
2996 
2997  /** \brief Const vector range function template for <tt>std::vector</tt>
2998 
2999  The element with index \c start in the original vector
3000  will become the first argument in the new vector, and
3001  the new vector will have size <tt>last-start</tt> .
3002 
3003  \note In this case, the return type is the same as the
3004  type of the first parameter.
3005  \note Unlike the ublas and pointer cases, this forces
3006  a copy.
3007  */
3008  template<class dat_t> const std::vector<dat_t>
3009  vector_range_copy(const std::vector<dat_t> &v, size_t start,
3010  size_t last) {
3011  return std::vector<dat_t> (v.begin()+start,v.begin()+last);
3012  }
3013 
3014  //@}
3015 
3016 }
3017 
3018 #if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)
3019 
3020 #if defined (O2SCL_ARMA) || defined (DOXYGEN)
3021 #include <armadillo>
3022 namespace o2scl {
3023 
3024  /// \name Armadillo specializations
3025  //@{
3026  /// Armadillo version of \ref matrix_max()
3027  double matrix_max(const arma::mat &data);
3028 
3029  /// Armadillo version of \ref matrix_min()
3030  double matrix_min(const arma::mat &data);
3031 
3032  /// Armadillo version of \ref matrix_row()
3033  template<> arma::subview_row<double>
3034  matrix_row<arma::mat,arma::subview_row<double> >
3035  (arma::mat &M, size_t row);
3036 
3037  /// Armadillo version of \ref matrix_column()
3038  template<> arma::subview_col<double>
3039  matrix_column<arma::mat,arma::subview_col<double> >
3040  (arma::mat &M, size_t column);
3041  //@}
3042 
3043 }
3044 
3045 #endif
3046 
3047 #if defined (O2SCL_EIGEN) || defined (DOXYGEN)
3048 #include <eigen3/Eigen/Dense>
3049 
3050 namespace o2scl {
3051 
3052  /// \name Eigen specializations
3053  //@{
3054  /// Eigen version of \ref matrix_max()
3055  double matrix_max(const Eigen::MatrixXd &data);
3056 
3057  /// Eigen version of \ref matrix_min()
3058  double matrix_min(const Eigen::MatrixXd &data);
3059 
3060  /// Eigen version of \ref matrix_row()
3061  template<> Eigen::MatrixXd::RowXpr
3062  matrix_row<Eigen::MatrixXd,Eigen::MatrixXd::RowXpr>
3063  (Eigen::MatrixXd &M, size_t row);
3064 
3065  /// Eigen version of \ref matrix_column()
3066  template<> Eigen::MatrixXd::ColXpr
3067  matrix_column<Eigen::MatrixXd,Eigen::MatrixXd::ColXpr>
3068  (Eigen::MatrixXd &M, size_t column);
3069  //@}
3070 
3071 }
3072 
3073 #endif
3074 
3075 #else
3076 
3077 #include <o2scl/vector_special.h>
3078 
3079 // End of "#if defined (O2SCL_COND_FLAG) || defined (DOXYGEN)"
3080 #endif
3081 
3082 #ifdef DOXYGEN
3083 /** \brief Placeholder documentation of some related Boost objects
3084  */
3085 namespace boost {
3086  /** \brief Documentation of Boost::numeric objects
3087  */
3088  namespace numeric {
3089  /** \brief Documentation of uBlas objects
3090  */
3091  namespace ublas {
3092  /** \brief The default vector type from uBlas
3093 
3094  The uBlas types aren't documented here, but the full documentation
3095  is available at
3096  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3097 
3098  Internally in \o2, this is often typedef'd using
3099  \code
3100  typedef boost::numeric::ublas::vector<double> ubvector;
3101  typedef boost::numeric::ublas::vector<size_t> ubvector_size_t;
3102  typedef boost::numeric::ublas::vector<int> ubvector_int;
3103  \endcode
3104 
3105  This is documented in \ref vector.h .
3106  */
3107  template<class T, class A> class vector {
3108  };
3109  /** \brief The default matrix type from uBlas
3110 
3111  The uBlas types aren't documented here, but the full documentation
3112  is available at
3113  http://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm
3114 
3115  Internally in \o2, this is often typedef'd using
3116  \code
3117  typedef boost::numeric::ublas::matrix<double> ubmatrix;
3118  typedef boost::numeric::ublas::matrix<size_t> ubmatrix_size_t;
3119  typedef boost::numeric::ublas::matrix<int> ubmatrix_int;
3120  \endcode
3121 
3122  This is documented in \ref vector.h .
3123  */
3124  template<class T, class F, class A> class matrix {
3125  };
3126  }
3127  }
3128 }
3129 // End of "#ifdef DOXYGEN"
3130 #endif
3131 
3132 // End of "#ifndef O2SCL_VECTOR_H"
3133 #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:2763
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
void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest)
From a given vector, create a new vector by removing a specified element.
Definition: vector.h:2299
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
bool vector_is_finite(size_t n, vec_t &data)
Test if the first n elements of a vector are finite.
Definition: vector.h:2109
Placeholder documentation of some related Boost objects.
Definition: vector.h:3085
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:2856
Generic object which represents a column of a matrix.
Definition: vector.h:2500
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:2816
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:2444
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:2950
void vector_reverse(size_t n, vec_t &data)
Reverse the first n elements of a vector.
Definition: vector.h:2350
size_t last_
The end() iterator.
Definition: vector.h:2776
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:2331
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:2851
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:3124
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:2238
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:2137
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:2993
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:2886
matrix_row_gen(mat_t &m, size_t row)
Create a row object from row row of matrix m.
Definition: vector.h:2455
vec_t & v_
A reference to the original vector.
Definition: vector.h:2770
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:2904
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:2845
dat_t * vector_range(dat_t *v, size_t start, size_t last)
Vector range function for pointers.
Definition: vector.h:2593
#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:2868
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:2575
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:2793
mat_column_t matrix_column(mat_t &M, size_t column)
Construct a column of a matrix.
Definition: vector.h:2486
size_t start_
The index offset.
Definition: vector.h:2773
#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:2430
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:2528
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:2568
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:2194
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:2464
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:2827
const double & operator[](size_t i) const
Return a const reference ith element.
Definition: vector.h:2909
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:2165
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:2781
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:2759
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:2385
int vector_is_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are monotonic and increasing or decreasing.
Definition: vector.h:1976
int vector_is_strictly_monotonic(size_t n, vec_t &data)
Test if the first n elements of a vector are strictly monotonic and determine if they are increasing ...
Definition: vector.h:2058
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
void vector_set_all(size_t N, vec_t &src, data_t &val)
Set the first N entries in a vector to a particular value.
Definition: vector.h:2280
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:2848
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:2603
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:2459
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:3107
size_t size() const
Return the vector size.
Definition: vector.h:2811
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).