/usr/include/rheolef/ic0.h is in librheolef-dev 5.93-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 | #ifndef _SKIT_IC0_H
#define _SKIT_IC0_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
///
/// =========================================================================
/*Class:ic0
NAME: @code{ic0} - incomplete Choleski factorization
@clindex ssr
@clindex vec
@cindex sparse matrix
DESCRIPTION:
@noindent
The class implements the icomplete Choleski factorization IC0(A)
when $A$ is a square symmetric matrix stored in CSR format.
@example
csr<double> a = ...;
ic0<double> fact(a);
@end example
AUTHOR:
Pierre Saramito
| Pierre.Saramito@imag.fr
LMC-IMAG, 38041 Grenoble cedex 9, France
DATE: 12 november 2009
End:
*/
#include "rheolef/csr.h"
namespace rheolef {
//<ic0:
template <class T>
class basic_ic0 : public csr<T> {
public:
typedef typename csr<T>::size_type size_type;
typedef T element_type;
// constructors:
basic_ic0 ();
explicit basic_ic0 (const csr<T>& a);
// solver:
void inplace_solve (vec<T>& x) const;
void solve (const vec<T>& b, vec<T>& x) const;
vec<T> solve (const vec<T>& b) const;
// accessors:
size_type nrow () const { return csr<T>::nrow(); }
size_type ncol () const { return csr<T>::ncol(); }
size_type nnz () const { return csr<T>::nnz(); }
};
template <class T>
basic_ic0<T> ic0 (const csr<T>& a);
//>ic0:
}// namespace rheolef
// ------------------------------------------------------------------
// inlined implementations
// ------------------------------------------------------------------
# include "rheolef/ic0-algo.h"
namespace rheolef {
template<class T>
inline
basic_ic0<T>::basic_ic0 ()
: csr<T>()
{
}
template<class T>
inline
basic_ic0<T>::basic_ic0 (const csr<T>& a)
: csr<T>()
{
size_t nnz_extra_diag_upper = csr_nnz_extra_diag_upper (a.nrow(), a.ia().begin(), a.ja().begin());
csr<T>::resize (a.nrow(), a.ncol(), nnz_extra_diag_upper + a.nrow());
#ifdef TO_CLEAN
cerr << "A = " << matlab << a << endl;
#endif // TO_CLEAN
csr_extract_upper (a.nrow(),
a.ia().begin(), a.ja().begin(), a.a().begin(),
csr<T>::ia().begin(), csr<T>::ja().begin(), csr<T>::a().begin());
#ifdef TO_CLEAN
cerr << "CSC_A = " << matlab << csr<T>(*this) << endl;
#endif // TO_CLEAN
csc_ic0 (a.nrow(),
csr<T>::ia().begin(), csr<T>::ja().begin(), csr<T>::a().begin());
#ifdef TO_CLEAN
cerr << "IC0_A = " << matlab << csr<T>(*this) << endl;
#endif // TO_CLEAN
}
template <class T>
inline
basic_ic0<T> ic0 (const csr<T>& a)
{
return basic_ic0<T>(a);
}
template<class T>
inline
void
basic_ic0<T>::inplace_solve (vec<T>& x) const
{
fatal_macro ("not yet");
}
template<class T>
inline
void
basic_ic0<T>::solve (const vec<T>& b, vec<T>& x) const
{
std::copy (b.begin(), b.end(), x.begin());
// solve L*y = b (use x for y)
csc_lower_diag_solve (nrow(),
csr<T>::ia().begin(), csr<T>::ja().begin(), csr<T>::a().begin(),
x.begin());
// solve L*x = y (use x for y)
csr_upper_diag_solve (nrow(),
csr<T>::ia().begin(), csr<T>::ja().begin(), csr<T>::a().begin(),
x.begin(), x.begin());
}
template<class T>
inline
vec<T>
basic_ic0<T>::solve (const vec<T>& b) const
{
vec<T> x (b.size());
solve (b,x);
#ifdef TO_CLEAN
cerr << "b = " << matlab << b << endl;
cerr << "x = " << matlab << x << endl;
#endif // TO_CLEAN
return x;
}
}// namespace rheolef
#endif // _SKIT_IC0_H
|