lu.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file lu.h
24  \brief Header wrapper for \ref lu_base.h
25 */
26 #ifndef O2SCL_LU_H
27 #define O2SCL_LU_H
28 
29 #include <o2scl/err_hnd.h>
30 #include <o2scl/permutation.h>
31 #include <o2scl/cblas.h>
32 
33 namespace o2scl_linalg {
34 
35 #define O2SCL_IX(V,i) V[i]
36 #define O2SCL_IX2(M,i,j) M(i,j)
37 #include <o2scl/lu_base.h>
38 #undef O2SCL_IX
39 #undef O2SCL_IX2
40 
41 }
42 
43 namespace o2scl_linalg_bracket {
44 
45  /** \brief Specialized version of LU_decomp for C-style 2D arrays
46 
47  \note Note that \c N and \c n must be equal, by no checking
48  is done to ensure that this is the case
49  */
50  template<size_t N>
51  int LU_decomp_array_2d(const size_t n, double A[][N],
52  o2scl::permutation &p, int &signum) {
53 
54  size_t i, j, k;
55 
56  signum=1;
57  p.init();
58 
59  for (j = 0; j < N - 1; j++) {
60 
61  /* Find maximum in the j-th column */
62  double ajj, max = fabs(A[j][j]);
63  size_t i_pivot = j;
64 
65  for (i = j + 1; i < N; i++) {
66  double aij = fabs (A[i][j]);
67 
68  if (aij > max) {
69  max = aij;
70  i_pivot = i;
71  }
72  }
73 
74  if (i_pivot != j) {
75 
76  // Swap rows j and i_pivot
77  double temp;
78  double *r1=&(A[j][0]);
79  double *r2=&(A[i_pivot][0]);
80  for (k=0;k<N;k++) {
81  temp=r1[k];
82  r1[k]=r2[k];
83  r2[k]=temp;
84  }
85  p.swap(j,i_pivot);
86  signum=-signum;
87  }
88 
89  ajj = A[j][j];
90 
91  if (ajj != 0.0) {
92  for (i = j + 1; i < N; i++) {
93  double aij = A[i][j] / ajj;
94  A[i][j]=aij;
95  for (k = j + 1; k < N; k++) {
96  double aik = A[i][k];
97  double ajk = A[j][k];
98  A[i][k]=aik - aij * ajk;
99  }
100  }
101  }
102  }
103 
104  return o2scl::success;
105  }
106 
107 #define O2SCL_IX(V,i) V[i]
108 #define O2SCL_IX2(M,i,j) M[i][j]
109 #include <o2scl/lu_base.h>
110 #undef O2SCL_IX
111 #undef O2SCL_IX2
112 
113 }
114 
115 #endif
int LU_decomp_array_2d(const size_t n, double A[][N], o2scl::permutation &p, int &signum)
Specialized version of LU_decomp for C-style 2D arrays.
Definition: lu.h:51
A class for representing permutations.
Definition: permutation.h:70
The namespace for linear algebra classes and functions with operator()
Definition: bidiag.h:46
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
int init()
Initialize permutation to the identity.
Definition: permutation.h:203
int swap(const size_t i, const size_t j)
Swap two elements of a permutation.
Definition: permutation.h:249
Success.
Definition: err_hnd.h:47

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