// forward substitution
for (int i=0; i < n; i++)
{
- Real sum (0.0);
- for (int j=0; j < i; j++)
- sum += y (j) * L(i,j);
- y (i) = (rhs (i) - sum)/L(i,i);
+ Real sum (0.0);
+ for (int j=0; j < i; j++)
+ sum += y (j) * L(i,j);
+ y (i) = (rhs (i) - sum)/L(i,i);
}
for (int i=0; i < n; i++)
- y (i) /= D(i);
+ y (i) /= D(i);
// backward subst
Vector &x (out); // using input as return val.
for (int i=n-1; i >= 0; i--)
{
- Real sum (0.0);
- for (int j=i+1; j < n; j++)
- sum += L(j,i)*x (j);
- x (i) = (y (i) - sum)/L(i,i);
+ Real sum (0.0);
+ for (int j=i+1; j < n; j++)
+ sum += L(j,i)*x (j);
+ x (i) = (y (i) - sum)/L(i,i);
}
}
int b = L.band_i();
assert (n == L.dim());
- out.set_dim (n);
+ out.set_dim (n);
Vector y;
y.set_dim (n);
// forward substitution
for (int i=0; i < n; i++)
{
- Real sum (0.0);
- for (int j= 0 >? i - b; j < i; j++)
- sum += y (j) * L(i,j);
- y (i) = (rhs (i) - sum)/L(i,i);
+ Real sum (0.0);
+ for (int j= 0 >? i - b; j < i; j++)
+ sum += y (j) * L(i,j);
+ y (i) = (rhs (i) - sum)/L(i,i);
}
for (int i=0; i < n; i++)
- y (i) /= D(i);
+ y (i) /= D(i);
// backward subst
Vector &x (out); // using input as return val.
for (int i=n-1; i >= 0; i--)
{
- Real sum (0.0);
- for (int j=i+1; j <= i + b&&j < n ; j++)
- sum += L(j,i)*x (j);
- x (i) = (y (i) - sum)/L(i,i);
+ Real sum (0.0);
+ for (int j=i+1; j <= i + b&&j < n ; j++)
+ sum += L(j,i)*x (j);
+ x (i) = (y (i) - sum)/L(i,i);
}
}
{
if (L.band_b())
{
- band_matrix_solve (x,rhs);
+ band_matrix_solve (x,rhs);
}
else
- full_matrix_solve (x,rhs);
+ full_matrix_solve (x,rhs);
}
Vector
Choleski_decomposition::full_matrix_decompose (Matrix const & P)
{
- int n = P.dim();
+ int n = P.dim();
L.unit();
for (int k= 0; k < n; k++)
{
- for (int j = 0; j < k; j++)
- {
- Real sum (0.0);
- for (int l=0; l < j; l++)
- sum += L(k,l)*L(j,l)*D(l);
- L(k,j) = (P(k,j) - sum)/D(j);
- }
- Real sum=0.0;
+ for (int j = 0; j < k; j++)
+ {
+ Real sum (0.0);
+ for (int l=0; l < j; l++)
+ sum += L(k,l)*L(j,l)*D(l);
+ L(k,j) = (P(k,j) - sum)/D(j);
+ }
+ Real sum=0.0;
- for (int l=0; l < k; l++)
- sum += sqr (L(k,l))*D(l);
- Real d = P(k,k) - sum;
- D(k) = d;
+ for (int l=0; l < k; l++)
+ sum += sqr (L(k,l))*D(l);
+ Real d = P(k,k) - sum;
+ D(k) = d;
}
}
for (int i= 0; i < n; i++)
{
- for (int j = 0 >? i - b; j < i; j++)
- {
- Real sum (0.0);
- for (int l=0 >? i - b; l < j; l++)
- sum += L(i,l)*L(j,l)*D(l);
- L(i,j) = (P(i,j) - sum)/D(j);
- }
- Real sum=0.0;
+ for (int j = 0 >? i - b; j < i; j++)
+ {
+ Real sum (0.0);
+ for (int l=0 >? i - b; l < j; l++)
+ sum += L(i,l)*L(j,l)*D(l);
+ L(i,j) = (P(i,j) - sum)/D(j);
+ }
+ Real sum=0.0;
- for (int l=0 >? i - b; l < i; l++)
- sum += sqr (L(i,l))*D(l);
- Real d = P(i,i) - sum;
- D(i) = d;
+ for (int l=0 >? i - b; l < i; l++)
+ sum += sqr (L(i,l))*D(l);
+ Real d = P(i,i) - sum;
+ D(i) = d;
}
L.try_set_band();
assert (L.band_i() == P.band_i ());
*/
Choleski_decomposition::Choleski_decomposition (Matrix const & P)
- : L(P.dim()), D(P.dim ())
+ : L(P.dim()), D(P.dim ())
{
#ifdef PARANOID
assert ((P-P.transposed()).norm ()/P.norm () < EPS);
#endif
if (P.band_b())
- band_matrix_decompose (P);
+ band_matrix_decompose (P);
else
- full_matrix_decompose (P);
+ full_matrix_decompose (P);
#ifdef PARANOID
Vector inv (n);
for (int i = 0; i < n; i++)
{
- e_i.set_unit (i);
- solve (inv, e_i);
- for (int j = 0 ; j<n; j++)
- invm (i,j) = inv (j);
+ e_i.set_unit (i);
+ solve (inv, e_i);
+ for (int j = 0 ; j<n; j++)
+ invm (i,j) = inv (j);
}
#ifdef PARANOID
Directed_graph_node::copy_edges_out (Directed_graph_node const &s)
{
for (int i=0; i < s.edge_out_l_arr_.size(); i++)
- add (s.edge_out_l_arr_[i]);
+ add (s.edge_out_l_arr_[i]);
}
void
#ifndef NDEBUG
for (int i=0; i < edge_out_l_arr_.size(); i++)
{
- assert (edge_out_l_arr_[i]->
- edge_in_l_arr_.find_l (this));
+ assert (edge_out_l_arr_[i]->
+ edge_in_l_arr_.find_l (this));
}
for (int i=0; i < edge_in_l_arr_.size(); i++)
- assert (edge_in_l_arr_[i]->contains_b (this));
+ assert (edge_in_l_arr_[i]->contains_b (this));
#endif
}
PARANOID_OK();
for (int i=0; i < edge_out_l_arr_.size();)
{
- if (edge_out_l_arr_[i]== d_l)
- remove_edge_out_idx (i);
- else
- i++;
+ if (edge_out_l_arr_[i]== d_l)
+ remove_edge_out_idx (i);
+ else
+ i++;
}
PARANOID_OK();
}
Directed_graph_node::unlink()
{
#ifdef PARANOID
- PARANOID_OK();
+ PARANOID_OK();
- Link_array<Directed_graph_node> t = edge_out_l_arr_;
- t.concat (edge_in_l_arr_);
+ Link_array<Directed_graph_node> t = edge_out_l_arr_;
+ t.concat (edge_in_l_arr_);
#endif
- while (edge_out_l_arr_.size())
- remove_edge_out_idx (0);
+ while (edge_out_l_arr_.size())
+ remove_edge_out_idx (0);
- while (edge_in_l_arr_.size())
- remove_edge_in (edge_in_l_arr_[0]);
+ while (edge_in_l_arr_.size())
+ remove_edge_in (edge_in_l_arr_[0]);
#ifdef PARANOID
- for (int i =0; i < t.size(); i++)
- t[i]->OK();
+ for (int i =0; i < t.size(); i++)
+ t[i]->OK();
#endif
}
{
PARANOID_OK();
if (!dep_l)
- return ;
+ return ;
dep_l->edge_in_l_arr_.push (this);
edge_out_l_arr_.push (dep_l);
PARANOID_OK();
return -1;
char const* me = strh_.ch_C();
- char const* p = memrchr (me, length_i(), c);
+ char const* p = memrchr ((Byte*)me, length_i(), c);
if (p)
return p - me;
return -1;
}
/**
- find the substring.
+ find a substring.
@return
1 index of leftmost occurrence of #searchfor#
String::index_i (String searchfor) const
{
char const* me = strh_.ch_C();
- char const* p = (char const *) memmem (
- me, length_i(), searchfor.ch_C(), searchfor.length_i ());
+
+ char const* p = (char const *)
+ memmem (me, length_i(), searchfor.ch_C(), searchfor.length_i ());
if (p)
return p - me;
return String_convert::dec2_f (*this);
}
-