This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch master
in repository getfem.
The following commit(s) were added to refs/heads/master by this push:
new 956a40ef Remove leftovers and redundant includes, fix whitespace and
other minor changes
956a40ef is described below
commit 956a40eff0cb8a34dba51cce483ccd5d6f27e4c7
Author: Konstantinos Poulios <logar...@gmail.com>
AuthorDate: Tue Dec 26 18:35:16 2023 +0100
Remove leftovers and redundant includes, fix whitespace and other minor
changes
---
configure.ac | 8 +-
contrib/continuum_mechanics/Makefile.am | 2 -
interface/src/gf_spmat_get.cc | 2 +-
src/Makefile.am | 2 +-
src/dal_singleton.cc | 6 +-
src/getfem_superlu.cc | 2 +-
src/gmm/gmm_MUMPS_interface.h | 2 +-
src/gmm/gmm_dense_qr.h | 11 +-
src/gmm/gmm_solver_Schwarz_additive.h | 457 ++++++++++++++++----------------
9 files changed, 244 insertions(+), 248 deletions(-)
diff --git a/configure.ac b/configure.ac
index 018bfd86..4fcf6c1f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -58,7 +58,7 @@ dnl set the optimization level
dnl --------------------------
AC_ARG_WITH(optimization,
- AC_HELP_STRING([--with-optimization=FLAG],[Set the optimization
level (-O3 by default)]),
+ AS_HELP_STRING([--with-optimization=FLAG],[Set the optimization
level (-O3 by default)]),
[with_optimization=$withval],
[with_optimization='-O3']
)
@@ -647,12 +647,12 @@ LIBS="$acx_blas_save_LIBS"
# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
if test x"$acx_blas_ok" = xyes; then
- echo "OK, You have working BLAS libs ! Using $BLAS_LIBS" ;
HAVE_VENDOR_BLAS=1
+ echo "OK, You have working BLAS libs ! Using $BLAS_LIBS" ; HAVE_VENDOR_BLAS=1
else
- echo " *** YOU DONT HAVE BLAS! *** Using a cheap replacement" ;
HAVE_VENDOR_BLAS=0
+ echo " *** YOU DONT HAVE BLAS! *** Using a cheap replacement" ;
HAVE_VENDOR_BLAS=0
fi
-dnl ACX_BLAS([ echo "OK, You have working BLAS libs !"; HAVE_VENDOR_BLAS=1 ],
[echo "YOU DONT HAVE BLAS! Using a cheap replacement" ; HAVE_VENDOR_BLAS=0])
+dnl ACX_BLAS([ echo "OK, You have working BLAS libs !"; HAVE_VENDOR_BLAS=1
],[echo "YOU DONT HAVE BLAS! Using a cheap replacement" ; HAVE_VENDOR_BLAS=0])
LIBS="$LIBS $BLAS_LIBS"
CPPFLAGS="$CPPFLAGS -DGMM_USES_BLAS"
diff --git a/contrib/continuum_mechanics/Makefile.am
b/contrib/continuum_mechanics/Makefile.am
index dc5d9fc1..879ec8a9 100644
--- a/contrib/continuum_mechanics/Makefile.am
+++ b/contrib/continuum_mechanics/Makefile.am
@@ -22,8 +22,6 @@ EXTRA_DIST = \
check_PROGRAMS =
-CLEANFILES =
-
if BUILDPYTHON
TESTS = plasticity_fin_strain_lin_hardening_plane_strain.py
diff --git a/interface/src/gf_spmat_get.cc b/interface/src/gf_spmat_get.cc
index ad59332a..44d75228 100644
--- a/interface/src/gf_spmat_get.cc
+++ b/interface/src/gf_spmat_get.cc
@@ -403,7 +403,7 @@ void gf_spmat_get(getfemint::mexargs_in& m_in,
);
-#if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H)
+#if defined(GMM_USES_MUMPS)
/*@GET @CELL{mantissa_r, mantissa_i, exponent} = ('determinant')
returns the matrix determinant calculated using MUMPS.@*/
sub_command
diff --git a/src/Makefile.am b/src/Makefile.am
index 16d06fcc..f90d9dbc 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -256,6 +256,6 @@ libgetfem_la_LDFLAGS = ${LIBTOOL_VERSION_INFO}
libgetfem_la_LIBADD = @SUPERLU_LIBS@ @MUMPS_LIBS@
AM_CPPFLAGS = -I$(top_srcdir)/src -I../src -I$(top_srcdir)
-CLEANFILES = ii_files/* *.o.d getfem/getfem_im_list.h
+CLEANFILES = ii_files/* *.o.d getfem/getfem_im_list.h
DISTCLEANFILES = getfem/getfem_arch_config.h gmm/gmm_arch_config.h
diff --git a/src/dal_singleton.cc b/src/dal_singleton.cc
index 183f093c..91a87825 100644
--- a/src/dal_singleton.cc
+++ b/src/dal_singleton.cc
@@ -20,10 +20,6 @@
===========================================================================*/
#include "getfem/dal_singleton.h"
-#include <algorithm>
-#include "gmm/gmm.h"
-#include "getfem/getfem_omp.h"
-
namespace dal {
@@ -64,4 +60,4 @@ namespace dal {
}
}
-}/* end of namespace dal
*/
\ No newline at end of file
+}/* end of namespace dal
*/
diff --git a/src/getfem_superlu.cc b/src/getfem_superlu.cc
index 4e7c3538..fe7d2803 100644
--- a/src/getfem_superlu.cc
+++ b/src/getfem_superlu.cc
@@ -23,7 +23,7 @@
typedef int int_t;
-/* because SRC/util.h defines TRUE and FALSE ... */
+/* because slu_util.h defines TRUE and FALSE ... */
#ifdef TRUE
# undef TRUE
#endif
diff --git a/src/gmm/gmm_MUMPS_interface.h b/src/gmm/gmm_MUMPS_interface.h
index 968f18d2..f912c145 100644
--- a/src/gmm/gmm_MUMPS_interface.h
+++ b/src/gmm/gmm_MUMPS_interface.h
@@ -35,7 +35,7 @@
@date December 8, 2005.
@brief Interface with MUMPS (LU direct solver for sparse matrices).
*/
-#if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H)
+#if defined(GMM_USES_MUMPS)
#ifndef GMM_MUMPS_INTERFACE_H
#define GMM_MUMPS_INTERFACE_H
diff --git a/src/gmm/gmm_dense_qr.h b/src/gmm/gmm_dense_qr.h
index 4bdde208..e046a155 100644
--- a/src/gmm/gmm_dense_qr.h
+++ b/src/gmm/gmm_dense_qr.h
@@ -573,14 +573,15 @@ namespace gmm {
typedef typename linalg_traits<MAT1>::value_type T;
typedef typename number_traits<T>::magnitude_type R;
- size_type n(mat_nrows(A)), q(0), q_old, p(0), ite(0), its(0);
+ size_type n(mat_nrows(A)), q(0), p(0);
dense_matrix<T> H(n,n);
- sub_interval SUBK(0,0);
-
gmm::copy(A, H);
+
Hessenberg_reduction(H, Q, compvect);
qr_stop_criterion(H, p, q, tol);
+ size_type ite(0), its(0);
+ sub_interval SUBK(0,0);
while (q < n) {
sub_interval SUBI(p, n-p-q), SUBJ(0, mat_ncols(Q));
if (compvect) SUBK = SUBI;
@@ -589,7 +590,7 @@ namespace gmm {
Wilkinson_double_shift_qr_step(sub_matrix(H, SUBI),
sub_matrix(Q, SUBJ, SUBK),
tol, (its == 10 || its == 20), compvect);
- q_old = q;
+ size_type q_old(q);
qr_stop_criterion(H, p, q, tol*R(2));
if (q != q_old) its = 0;
++its; ++ite;
@@ -723,7 +724,7 @@ namespace gmm {
typedef typename linalg_traits<MAT1>::value_type T;
typedef typename number_traits<T>::magnitude_type R;
- size_type n = mat_nrows(A), q = 0, p, ite = 0;
+ size_type n(mat_nrows(A)), q(0), p(0), ite(0);
dense_matrix<T> Tri(n, n);
gmm::copy(A, Tri);
diff --git a/src/gmm/gmm_solver_Schwarz_additive.h
b/src/gmm/gmm_solver_Schwarz_additive.h
index 5585e027..a37bbef4 100644
--- a/src/gmm/gmm_solver_Schwarz_additive.h
+++ b/src/gmm/gmm_solver_Schwarz_additive.h
@@ -36,7 +36,7 @@
*/
#ifndef GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
-#define GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
+#define GMM_SOLVERS_SCHWARZ_ADDITIVE_H__
#include "gmm_kernel.h"
#include "gmm_superlu_interface.h"
@@ -46,9 +46,9 @@
#include "gmm_solver_qmr.h"
namespace gmm {
-
+
/* ******************************************************************** */
- /* Additive Schwarz interfaced local solvers */
+ /* Additive Schwarz interfaced local solvers */
/* ******************************************************************** */
struct using_cg {};
@@ -62,24 +62,24 @@ namespace gmm {
static const APrecond &transform(const P &PP) { return PP; }
};
- template <typename Matrix1, typename Precond, typename Vector>
+ template <typename Matrix1, typename Precond, typename Vector>
void AS_local_solve(using_cg, const Matrix1 &A, Vector &x, const Vector &b,
- const Precond &P, iteration &iter)
+ const Precond &P, iteration &iter)
{ cg(A, x, b, P, iter); }
- template <typename Matrix1, typename Precond, typename Vector>
+ template <typename Matrix1, typename Precond, typename Vector>
void AS_local_solve(using_gmres, const Matrix1 &A, Vector &x,
- const Vector &b, const Precond &P, iteration &iter)
+ const Vector &b, const Precond &P, iteration &iter)
{ gmres(A, x, b, P, 100, iter); }
-
- template <typename Matrix1, typename Precond, typename Vector>
+
+ template <typename Matrix1, typename Precond, typename Vector>
void AS_local_solve(using_bicgstab, const Matrix1 &A, Vector &x,
- const Vector &b, const Precond &P, iteration &iter)
+ const Vector &b, const Precond &P, iteration &iter)
{ bicgstab(A, x, b, P, iter); }
- template <typename Matrix1, typename Precond, typename Vector>
+ template <typename Matrix1, typename Precond, typename Vector>
void AS_local_solve(using_qmr, const Matrix1 &A, Vector &x,
- const Vector &b, const Precond &P, iteration &iter)
+ const Vector &b, const Precond &P, iteration &iter)
{ qmr(A, x, b, P, iter); }
#if defined(GMM_USES_SUPERLU)
@@ -94,18 +94,18 @@ namespace gmm {
static const APrecond &transform(const APrecond &PP) { return PP; }
};
- template <typename Matrix1, typename Precond, typename Vector>
+ template <typename Matrix1, typename Precond, typename Vector>
void AS_local_solve(using_superlu, const Matrix1 &, Vector &x,
- const Vector &b, const Precond &P, iteration &iter)
+ const Vector &b, const Precond &P, iteration &iter)
{ P.solve(x, b); iter.set_iteration(1); }
#endif
/* ******************************************************************** */
- /* Additive Schwarz Linear system */
+ /* Additive Schwarz Linear system */
/* ******************************************************************** */
template <typename Matrix1, typename Matrix2, typename Precond,
- typename local_solver>
+ typename local_solver>
struct add_schwarz_mat{
typedef typename linalg_traits<Matrix1>::value_type value_type;
@@ -117,34 +117,34 @@ namespace gmm {
mutable size_type itebilan;
mutable std::vector<std::vector<value_type> > gi, fi;
std::vector<typename actual_precond<Precond, local_solver,
- Matrix1>::APrecond> precond1;
+ Matrix1>::APrecond> precond1;
void init(const Matrix1 &A_, const std::vector<Matrix2> &vB_,
- iteration iter_, const Precond &P, double residual_);
+ iteration iter_, const Precond &P, double residual_);
- add_schwarz_mat(void) {}
+ add_schwarz_mat() {}
add_schwarz_mat(const Matrix1 &A_, const std::vector<Matrix2> &vB_,
- iteration iter_, const Precond &P, double residual_)
+ iteration iter_, const Precond &P, double residual_)
{ init(A_, vB_, iter_, P, residual_); }
};
template <typename Matrix1, typename Matrix2, typename Precond,
- typename local_solver>
+ typename local_solver>
void add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver>::init(
const Matrix1 &A_, const std::vector<Matrix2> &vB_,
iteration iter_, const Precond &P, double residual_) {
vB = &vB_; A = &A_; iter = iter_;
residual = residual_;
-
+
size_type nb_sub = vB->size();
vAloc.resize(nb_sub);
gi.resize(nb_sub); fi.resize(nb_sub);
precond1.resize(nb_sub);
std::fill(precond1.begin(), precond1.end(),
- actual_precond<Precond, local_solver, Matrix1>::transform(P));
+ actual_precond<Precond, local_solver, Matrix1>::transform(P));
itebilan = 0;
-
+
if (iter.get_noisy()) cout << "Init pour sub dom ";
#ifdef GMM_USES_MPI
int size,tranche,borne_sup,borne_inf,rank,tag1=11,tag2=12,tag3=13,sizepr =
0;
@@ -167,47 +167,47 @@ namespace gmm {
int next = (rank + 1) % size;
int previous = (rank + size - 1) % size;
//communication of local information on ring pattern
- //Each process receive Nproc-1 contributions
+ //Each process receive Nproc-1 contributions
for (int nproc = 0; nproc < size; ++nproc) {
- for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i)
{
-// for (size_type i = 0; i < nb_sub/size; ++i) {
-// for (size_type i = 0; i < nb_sub; ++i) {
- // size_type i=(rank+size*(j-1)+nb_sub)%nb_sub;
+ for (size_type i = size_type(borne_inf); i < size_type(borne_sup); ++i) {
+// for (size_type i = 0; i < nb_sub/size; ++i) {
+// for (size_type i = 0; i < nb_sub; ++i) {
+ // size_type i=(rank+size*(j-1)+nb_sub)%nb_sub;
- cout << "Sous domaines " << i << " : " << mat_ncols((*vB)[i]) << endl;
+ cout << "Sous domaines " << i << " : " << mat_ncols((*vB)[i]) << endl;
#else
- for (size_type i = 0; i < nb_sub; ++i) {
+ for (size_type i = 0; i < nb_sub; ++i) {
#endif
-
- if (iter.get_noisy()) cout << i << " " << std::flush;
- Matrix2 Maux(mat_ncols((*vB)[i]), mat_nrows((*vB)[i]));
-
+ if (iter.get_noisy())
+ cout << i << " " << std::flush;
+ Matrix2 Maux(mat_ncols((*vB)[i]), mat_nrows((*vB)[i]));
+
#ifdef GMM_USES_MPI
- Matrix2 Maux2(mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
- if (nproc == 0) {
- gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
- gmm::clear(vAloc[i]);
- }
- gmm::mult(gmm::transposed((*vB)[i]), Acsr, Maux);
- gmm::mult(Maux, (*vB)[i], Maux2);
- gmm::add(Maux2, vAloc[i]);
+ Matrix2 Maux2(mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
+ if (nproc == 0) {
+ gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
+ gmm::clear(vAloc[i]);
+ }
+ gmm::mult(gmm::transposed((*vB)[i]), Acsr, Maux);
+ gmm::mult(Maux, (*vB)[i], Maux2);
+ gmm::add(Maux2, vAloc[i]);
#else
- gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
- gmm::mult(gmm::transposed((*vB)[i]), *A, Maux);
- gmm::mult(Maux, (*vB)[i], vAloc[i]);
+ gmm::resize(vAloc[i], mat_ncols((*vB)[i]), mat_ncols((*vB)[i]));
+ gmm::mult(gmm::transposed((*vB)[i]), *A, Maux);
+ gmm::mult(Maux, (*vB)[i], vAloc[i]);
#endif
#ifdef GMM_USES_MPI
- if (nproc == size - 1 ) {
+ if (nproc == size - 1 ) {
#endif
- precond1[i].build_with(vAloc[i]);
- gmm::resize(fi[i], mat_ncols((*vB)[i]));
- gmm::resize(gi[i], mat_ncols((*vB)[i]));
+ precond1[i].build_with(vAloc[i]);
+ gmm::resize(fi[i], mat_ncols((*vB)[i]));
+ gmm::resize(gi[i], mat_ncols((*vB)[i]));
#ifdef GMM_USES_MPI
- }
+ }
#else
- }
+ }
#endif
#ifdef GMM_USES_MPI
}
@@ -223,8 +223,8 @@ namespace gmm {
MPI_Sendrecv(&(Acsr.ir[0]), Acsr.jc[sizeA], MPI_INT, next, tag1,
&(Acsrtemp.ir[0]), Acsrtemp.jc[sizeA], MPI_INT, previous,
tag1,
MPI_COMM_WORLD, &status);
-
- MPI_Sendrecv(&(Acsr.pr[0]), Acsr.jc[sizeA], mpi_type(value_type()),
next, tag3,
+
+ MPI_Sendrecv(&(Acsr.pr[0]), Acsr.jc[sizeA], mpi_type(value_type()),
next, tag3,
&(Acsrtemp.pr[0]), Acsrtemp.jc[sizeA],
mpi_type(value_type()), previous, tag3,
MPI_COMM_WORLD, &status);
gmm::copy(Acsrtemp, Acsr);
@@ -235,11 +235,11 @@ namespace gmm {
#endif
if (iter.get_noisy()) cout << "\n";
}
-
+
template <typename Matrix1, typename Matrix2, typename Precond,
- typename Vector2, typename Vector3, typename local_solver>
+ typename Vector2, typename Vector3, typename local_solver>
void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
- const Vector2 &p, Vector3 &q) {
+ const Vector2 &p, Vector3 &q) {
size_type itebilan = 0;
#ifdef GMM_USES_MPI
static double tmult_tot = 0.0;
@@ -276,12 +276,12 @@ namespace gmm {
#endif
{
#ifdef GMM_USES_MPI
- // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
+ // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
#endif
- gmm::mult(gmm::transposed((*(M.vB))[i]), q, M.fi[i]);
+ gmm::mult(gmm::transposed((*(M.vB))[i]), q, M.fi[i]);
M.iter.init();
AS_local_solve(local_solver(), (M.vAloc)[i], (M.gi)[i],
- (M.fi)[i],(M.precond1)[i],M.iter);
+ (M.fi)[i],(M.precond1)[i],M.iter);
itebilan = std::max(itebilan, M.iter.get_iteration());
}
@@ -298,17 +298,17 @@ namespace gmm {
#else
for (size_type i = 0; i < M.gi.size(); ++i)
#endif
- {
+ {
#ifdef GMM_USES_MPI
- // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
-// gmm::mult((*(M.vB))[i], M.gi[i], qbis,qbis);
- gmm::mult((*(M.vB))[i], M.gi[i], qter);
- add(qter,qbis,qbis);
+ // size_type i=j; // (rank+size*(j-1)+nb_sub)%nb_sub;
+// gmm::mult((*(M.vB))[i], M.gi[i], qbis,qbis);
+ gmm::mult((*(M.vB))[i], M.gi[i], qter);
+ add(qter,qbis,qbis);
#else
- gmm::mult((*(M.vB))[i], M.gi[i], q, q);
+ gmm::mult((*(M.vB))[i], M.gi[i], q, q);
#endif
- }
+ }
#ifdef GMM_USES_MPI
//WARNING this add only if you use the ring pattern below
// need to do this below if using a n explicit ring pattern communication
@@ -326,25 +326,25 @@ namespace gmm {
// int next = (rank + 1) % size;
// int previous = (rank + size - 1) % size;
//communication of local information on ring pattern
- //Each process receive Nproc-1 contributions
+ //Each process receive Nproc-1 contributions
// if (size > 1) {
-// for (int nproc = 0; nproc < size-1; ++nproc)
+// for (int nproc = 0; nproc < size-1; ++nproc)
// {
-// MPI_Sendrecv(&(qbis[0]), gmm::vect_size(q), MPI_DOUBLE, next, tag1,
-// &(qter[0]), gmm::vect_size(q),MPI_DOUBLE,previous,tag1,
-// MPI_COMM_WORLD,&status);
-// gmm::copy(qter, qbis);
-// add(qbis,q,q);
+// MPI_Sendrecv(&(qbis[0]), gmm::vect_size(q), MPI_DOUBLE, next, tag1,
+// &(qter[0]), gmm::vect_size(q),MPI_DOUBLE,previous,tag1,
+// MPI_COMM_WORLD,&status);
+// gmm::copy(qter, qbis);
+// add(qbis,q,q);
// }
// }
MPI_Allreduce(&(qbis[0]), &(q[0]),gmm::vect_size(q), MPI_DOUBLE,
- MPI_SUM,MPI_COMM_WORLD);
+ MPI_SUM,MPI_COMM_WORLD);
t_final=MPI_Wtime();
t_tot += t_final-t_ref;
cout<<"["<< rank<<"] temps reduce Resol "<< t_final-t_ref << " t_tot = "
<< t_tot << endl;
-#endif
+#endif
if (M.iter.get_noisy() > 0) cout << "itebloc = " << itebilan << endl;
M.itebilan += itebilan;
@@ -352,61 +352,61 @@ namespace gmm {
}
template <typename Matrix1, typename Matrix2, typename Precond,
- typename Vector2, typename Vector3, typename local_solver>
+ typename Vector2, typename Vector3, typename local_solver>
void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
- const Vector2 &p, const Vector3 &q) {
+ const Vector2 &p, const Vector3 &q) {
mult(M, p, const_cast<Vector3 &>(q));
}
template <typename Matrix1, typename Matrix2, typename Precond,
- typename Vector2, typename Vector3, typename Vector4,
- typename local_solver>
+ typename Vector2, typename Vector3, typename Vector4,
+ typename local_solver>
void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
- const Vector2 &p, const Vector3 &p2, Vector4 &q)
+ const Vector2 &p, const Vector3 &p2, Vector4 &q)
{ mult(M, p, q); add(p2, q); }
template <typename Matrix1, typename Matrix2, typename Precond,
- typename Vector2, typename Vector3, typename Vector4,
- typename local_solver>
+ typename Vector2, typename Vector3, typename Vector4,
+ typename local_solver>
void mult(const add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &M,
- const Vector2 &p, const Vector3 &p2, const Vector4 &q)
+ const Vector2 &p, const Vector3 &p2, const Vector4 &q)
{ mult(M, p, const_cast<Vector4 &>(q)); add(p2, q); }
/* ******************************************************************** */
- /* Additive Schwarz interfaced global solvers */
+ /* Additive Schwarz interfaced global solvers */
/* ******************************************************************** */
template <typename ASM_type, typename Vect>
void AS_global_solve(using_cg, const ASM_type &ASM, Vect &x,
- const Vect &b, iteration &iter)
+ const Vect &b, iteration &iter)
{ cg(ASM, x, b, *(ASM.A), identity_matrix(), iter); }
template <typename ASM_type, typename Vect>
void AS_global_solve(using_gmres, const ASM_type &ASM, Vect &x,
- const Vect &b, iteration &iter)
+ const Vect &b, iteration &iter)
{ gmres(ASM, x, b, identity_matrix(), 100, iter); }
template <typename ASM_type, typename Vect>
void AS_global_solve(using_bicgstab, const ASM_type &ASM, Vect &x,
- const Vect &b, iteration &iter)
+ const Vect &b, iteration &iter)
{ bicgstab(ASM, x, b, identity_matrix(), iter); }
template <typename ASM_type, typename Vect>
void AS_global_solve(using_qmr,const ASM_type &ASM, Vect &x,
- const Vect &b, iteration &iter)
+ const Vect &b, iteration &iter)
{ qmr(ASM, x, b, identity_matrix(), iter); }
#if defined(GMM_USES_SUPERLU)
template <typename ASM_type, typename Vect>
void AS_global_solve(using_superlu, const ASM_type &, Vect &,
- const Vect &, iteration &) {
+ const Vect &, iteration &) {
GMM_ASSERT1(false, "You cannot use SuperLU as "
- "global solver in additive Schwarz meethod");
+ "global solver in additive Schwarz meethod");
}
#endif
-
+
/* ******************************************************************** */
- /* Linear Additive Schwarz method */
+ /* Linear Additive Schwarz method */
/* ******************************************************************** */
/* ref : Domain decomposition algorithms for the p-version finite */
/* element method for elliptic problems, Luca F. Pavarino, */
@@ -417,8 +417,8 @@ namespace gmm {
* with the same system.
*/
template <typename Matrix1, typename Matrix2,
- typename Vector2, typename Vector3, typename Precond,
- typename local_solver, typename global_solver>
+ typename Vector2, typename Vector3, typename Precond,
+ typename local_solver, typename global_solver>
void additive_schwarz(
add_schwarz_mat<Matrix1, Matrix2, Precond, local_solver> &ASM, Vector3 &u,
const Vector2 &f, iteration &iter, const global_solver&) {
@@ -453,11 +453,11 @@ namespace gmm {
gmm::mult(gmm::transposed((*(ASM.vB))[i]), f, ASM.fi[i]);
ASM.iter.init();
AS_local_solve(local_solver(), ASM.vAloc[i], ASM.gi[i], ASM.fi[i],
- ASM.precond1[i], ASM.iter);
+ ASM.precond1[i], ASM.iter);
ASM.itebilan = std::max(ASM.itebilan, ASM.iter.get_iteration());
#ifdef GMM_USES_MPI
gmm::mult((*(ASM.vB))[i], ASM.gi[i], gbis,gbis);
-#else
+#else
gmm::mult((*(ASM.vB))[i], ASM.gi[i], g, g);
#endif
}
@@ -466,7 +466,7 @@ namespace gmm {
double t_ref,t_final;
t_ref=MPI_Wtime();
MPI_Allreduce(&(gbis[0]), &(g[0]),gmm::vect_size(g), MPI_DOUBLE,
- MPI_SUM,MPI_COMM_WORLD);
+ MPI_SUM,MPI_COMM_WORLD);
t_final=MPI_Wtime();
cout<<"temps reduce init "<< t_final-t_ref<<endl;
#endif
@@ -487,13 +487,13 @@ namespace gmm {
* The ASM matrix represent the preconditionned linear system.
*/
template <typename Matrix1, typename Matrix2,
- typename Vector2, typename Vector3, typename Precond,
- typename local_solver, typename global_solver>
+ typename Vector2, typename Vector3, typename Precond,
+ typename local_solver, typename global_solver>
void additive_schwarz(const Matrix1 &A, Vector3 &u,
- const Vector2 &f, const Precond &P,
- const std::vector<Matrix2> &vB,
- iteration &iter, local_solver,
- global_solver) {
+ const Vector2 &f, const Precond &P,
+ const std::vector<Matrix2> &vB,
+ iteration &iter,
+ local_solver, global_solver) {
iter.set_rhsnorm(vect_norm2(f));
if (iter.get_rhsnorm() == 0.0) { gmm::clear(u); return; }
iteration iter2 = iter; iter2.reduce_noisy();
@@ -504,62 +504,62 @@ namespace gmm {
}
/* ******************************************************************** */
- /* Sequential Non-Linear Additive Schwarz method */
+ /* Sequential Non-Linear Additive Schwarz method */
/* ******************************************************************** */
/* ref : Nonlinearly Preconditionned Inexact Newton Algorithms, */
/* Xiao-Chuan Cai, David E. Keyes, */
- /* SIAM J. Sci. Comp. 24: p183-200. l */
+ /* SIAM J. Sci. Comp. 24: p183-200. l */
/* ******************************************************************** */
- template <typename Matrixt, typename MatrixBi>
+ template <typename Matrixt, typename MatrixBi>
class NewtonAS_struct {
-
+
public :
typedef Matrixt tangent_matrix_type;
typedef MatrixBi B_matrix_type;
typedef typename linalg_traits<Matrixt>::value_type value_type;
typedef std::vector<value_type> Vector;
-
- virtual size_type size(void) = 0;
+
+ virtual size_type size() = 0;
virtual const std::vector<MatrixBi> &get_vB() = 0;
-
+
virtual void compute_F(Vector &f, Vector &x) = 0;
virtual void compute_tangent_matrix(Matrixt &M, Vector &x) = 0;
// compute Bi^T grad(F(X)) Bi
virtual void compute_sub_tangent_matrix(Matrixt &Mloc, Vector &x,
- size_type i) = 0;
+ size_type i) = 0;
// compute Bi^T F(X)
virtual void compute_sub_F(Vector &fi, Vector &x, size_type i) = 0;
virtual ~NewtonAS_struct() {}
};
- template <typename Matrixt, typename MatrixBi>
+ template <typename Matrixt, typename MatrixBi>
struct AS_exact_gradient {
const std::vector<MatrixBi> &vB;
std::vector<Matrixt> vM;
std::vector<Matrixt> vMloc;
- void init(void) {
+ void init() {
for (size_type i = 0; i < vB.size(); ++i) {
- Matrixt aux(gmm::mat_ncols(vB[i]), gmm::mat_ncols(vM[i]));
- gmm::resize(vMloc[i], gmm::mat_ncols(vB[i]), gmm::mat_ncols(vB[i]));
- gmm::mult(gmm::transposed(vB[i]), vM[i], aux);
- gmm::mult(aux, vB[i], vMloc[i]);
+ Matrixt aux(gmm::mat_ncols(vB[i]), gmm::mat_ncols(vM[i]));
+ gmm::resize(vMloc[i], gmm::mat_ncols(vB[i]), gmm::mat_ncols(vB[i]));
+ gmm::mult(gmm::transposed(vB[i]), vM[i], aux);
+ gmm::mult(aux, vB[i], vMloc[i]);
}
}
AS_exact_gradient(const std::vector<MatrixBi> &vB_) : vB(vB_) {
vM.resize(vB.size()); vMloc.resize(vB.size());
for (size_type i = 0; i < vB.size(); ++i) {
- gmm::resize(vM[i], gmm::mat_nrows(vB[i]), gmm::mat_nrows(vB[i]));
+ gmm::resize(vM[i], gmm::mat_nrows(vB[i]), gmm::mat_nrows(vB[i]));
}
}
};
template <typename Matrixt, typename MatrixBi,
- typename Vector2, typename Vector3>
+ typename Vector2, typename Vector3>
void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
- const Vector2 &p, Vector3 &q) {
+ const Vector2 &p, Vector3 &q) {
gmm::clear(q);
typedef typename gmm::linalg_traits<Vector3>::value_type T;
std::vector<T> v(gmm::vect_size(p)), w, x;
@@ -577,26 +577,26 @@ namespace gmm {
}
template <typename Matrixt, typename MatrixBi,
- typename Vector2, typename Vector3>
+ typename Vector2, typename Vector3>
void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
- const Vector2 &p, const Vector3 &q) {
+ const Vector2 &p, const Vector3 &q) {
mult(M, p, const_cast<Vector3 &>(q));
}
template <typename Matrixt, typename MatrixBi,
- typename Vector2, typename Vector3, typename Vector4>
+ typename Vector2, typename Vector3, typename Vector4>
void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
- const Vector2 &p, const Vector3 &p2, Vector4 &q)
+ const Vector2 &p, const Vector3 &p2, Vector4 &q)
{ mult(M, p, q); add(p2, q); }
template <typename Matrixt, typename MatrixBi,
- typename Vector2, typename Vector3, typename Vector4>
+ typename Vector2, typename Vector3, typename Vector4>
void mult(const AS_exact_gradient<Matrixt, MatrixBi> &M,
- const Vector2 &p, const Vector3 &p2, const Vector4 &q)
+ const Vector2 &p, const Vector3 &p2, const Vector4 &q)
{ mult(M, p, const_cast<Vector4 &>(q)); add(p2, q); }
struct S_default_newton_line_search {
-
+
double conv_alpha, conv_r;
size_t it, itmax, glob_it;
@@ -607,9 +607,9 @@ namespace gmm {
double alpha_max_ratio_reached, r_max_ratio_reached;
size_type it_max_ratio_reached;
-
- double converged_value(void) { return conv_alpha; };
- double converged_residual(void) { return conv_r; };
+
+ double converged_value() { return conv_alpha; };
+ double converged_residual() { return conv_r; };
virtual void init_search(double r, size_t git, double = 0.0) {
alpha_min_ratio = 0.9;
@@ -623,7 +623,7 @@ namespace gmm {
count = 0;
max_ratio_reached = false;
}
- virtual double next_try(void) {
+ virtual double next_try() {
alpha_old = alpha;
if (alpha >= 0.4) alpha *= 0.5; else alpha *= alpha_mult; ++it;
return alpha_old;
@@ -631,45 +631,45 @@ namespace gmm {
virtual bool is_converged(double r, double = 0.0) {
// cout << "r = " << r << " alpha = " << alpha / alpha_mult << "
count_pat = " << count_pat << endl;
if (!max_ratio_reached && r < first_res * alpha_max_ratio) {
- alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
- it_max_ratio_reached = it; max_ratio_reached = true;
+ alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
+ it_max_ratio_reached = it; max_ratio_reached = true;
}
if (max_ratio_reached && r < r_max_ratio_reached * 0.5
- && r > first_res * 1.1 && it <= it_max_ratio_reached+1) {
- alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
- it_max_ratio_reached = it;
+ && r > first_res * 1.1 && it <= it_max_ratio_reached+1) {
+ alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r;
+ it_max_ratio_reached = it;
}
if (count == 0 || r < conv_r)
- { conv_r = r; conv_alpha = alpha_old; count = 1; }
+ { conv_r = r; conv_alpha = alpha_old; count = 1; }
if (conv_r < first_res) ++count;
if (r < first_res * alpha_min_ratio)
- { count_pat = 0; return true; }
+ { count_pat = 0; return true; }
if (count >= 5 || (alpha < alpha_min && max_ratio_reached)) {
- if (conv_r < first_res * 0.99) count_pat = 0;
- if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3)
- { conv_r=r_max_ratio_reached; conv_alpha=alpha_max_ratio_reached; }
- if (conv_r >= first_res * 0.9999) count_pat++;
- return true;
+ if (conv_r < first_res * 0.99) count_pat = 0;
+ if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3)
+ { conv_r=r_max_ratio_reached; conv_alpha=alpha_max_ratio_reached; }
+ if (conv_r >= first_res * 0.9999) count_pat++;
+ return true;
}
return false;
}
- S_default_newton_line_search(void) { count_pat = 0; }
+ S_default_newton_line_search() { count_pat = 0; }
};
-
+
template <typename Matrixt, typename MatrixBi, typename Vector,
- typename Precond, typename local_solver, typename global_solver>
+ typename Precond, typename local_solver, typename global_solver>
void Newton_additive_Schwarz(NewtonAS_struct<Matrixt, MatrixBi> &NS,
- const Vector &u_,
- iteration &iter, const Precond &P,
- local_solver, global_solver) {
+ const Vector &u_,
+ iteration &iter, const Precond &P,
+ local_solver, global_solver) {
Vector &u = const_cast<Vector &>(u_);
typedef typename linalg_traits<Vector>::value_type value_type;
typedef typename number_traits<value_type>::magnitude_type mtype;
typedef actual_precond<Precond, local_solver, Matrixt> chgt_precond;
-
+
double residual = iter.get_resmax();
S_default_newton_line_search internal_ls;
@@ -699,102 +699,103 @@ namespace gmm {
NS.compute_F(rhs, u);
mtype act_res=gmm::vect_norm2(rhs), act_res_new(0), precond_res = act_res;
mtype alpha;
-
+
while(!iter.finished(std::min(act_res, precond_res))) {
for (int SOR_step = 0; SOR_step >= 0; --SOR_step) {
- gmm::clear(rhs);
- for (size_type isd = 0; isd < NS.get_vB().size(); ++isd) {
- const MatrixBi &Bi = (NS.get_vB())[isd];
- size_type si = mat_ncols(Bi);
- gmm::resize(Mloc, si, si);
- xi.resize(si); xii.resize(si); fi.resize(si); di.resize(si);
-
- iternc.init();
- iternc.set_maxiter(30); // ?
- if (iternc.get_noisy())
- cout << "Non-linear local problem " << isd << endl;
- gmm::clear(xi);
- gmm::copy(u, x);
- NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
- mtype r = gmm::vect_norm2(fi), r_t(r);
- if (r > value_type(0)) {
- iternc.set_rhsnorm(std::max(r, mtype(1)));
- while(!iternc.finished(r)) {
- NS.compute_sub_tangent_matrix(Mloc, x, isd);
-
- PP.build_with(Mloc);
- iter3.init();
- AS_local_solve(local_solver(), Mloc, di, fi, PP, iter3);
-
- internal_ls.init_search(r, iternc.get_iteration());
- do {
- alpha = internal_ls.next_try();
- gmm::add(xi, gmm::scaled(di, -alpha), xii);
- gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
- NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
- r_t = gmm::vect_norm2(fi);
- } while (!internal_ls.is_converged(r_t));
-
- if (alpha != internal_ls.converged_value()) {
- alpha = internal_ls.converged_value();
- gmm::add(xi, gmm::scaled(di, -alpha), xii);
- gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
- NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
- r_t = gmm::vect_norm2(fi);
- }
- gmm::copy(x, vx[isd]); // for exact gradient
-
- if (iternc.get_noisy()) cout << "(step=" << alpha << ")\t";
- ++iternc; r = r_t; gmm::copy(xii, xi);
- }
- if (SOR_step) gmm::mult(Bi, gmm::scaled(xii, -1.0), u, u);
- gmm::mult(Bi, gmm::scaled(xii, -1.0), rhs, rhs);
- }
- }
- precond_res = gmm::vect_norm2(rhs);
- if (SOR_step) cout << "SOR step residual = " << precond_res << endl;
- if (precond_res < residual) break;
- cout << "Precond residual = " << precond_res << endl;
+ gmm::clear(rhs);
+ for (size_type isd = 0; isd < NS.get_vB().size(); ++isd) {
+ const MatrixBi &Bi = (NS.get_vB())[isd];
+ size_type si = mat_ncols(Bi);
+ gmm::resize(Mloc, si, si);
+ xi.resize(si); xii.resize(si); fi.resize(si); di.resize(si);
+
+ iternc.init();
+ iternc.set_maxiter(30); // ?
+ if (iternc.get_noisy())
+ cout << "Non-linear local problem " << isd << endl;
+ gmm::clear(xi);
+ gmm::copy(u, x);
+ NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
+ mtype r = gmm::vect_norm2(fi), r_t(r);
+ if (r > value_type(0)) {
+ iternc.set_rhsnorm(std::max(r, mtype(1)));
+ while(!iternc.finished(r)) {
+ NS.compute_sub_tangent_matrix(Mloc, x, isd);
+
+ PP.build_with(Mloc);
+ iter3.init();
+ AS_local_solve(local_solver(), Mloc, di, fi, PP, iter3);
+
+ internal_ls.init_search(r, iternc.get_iteration());
+ do {
+ alpha = internal_ls.next_try();
+ gmm::add(xi, gmm::scaled(di, -alpha), xii);
+ gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
+ NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
+ r_t = gmm::vect_norm2(fi);
+ } while (!internal_ls.is_converged(r_t));
+
+ if (alpha != internal_ls.converged_value()) {
+ alpha = internal_ls.converged_value();
+ gmm::add(xi, gmm::scaled(di, -alpha), xii);
+ gmm::mult(Bi, gmm::scaled(xii, -1.0), u, x);
+ NS.compute_sub_F(fi, x, isd); gmm::scale(fi, value_type(-1));
+ r_t = gmm::vect_norm2(fi);
+ }
+ gmm::copy(x, vx[isd]); // for exact gradient
+
+ if (iternc.get_noisy()) cout << "(step=" << alpha << ")\t";
+ ++iternc; r = r_t; gmm::copy(xii, xi);
+ }
+ if (SOR_step) gmm::mult(Bi, gmm::scaled(xii, -1.0), u, u);
+ gmm::mult(Bi, gmm::scaled(xii, -1.0), rhs, rhs);
+ }
+ }
+ precond_res = gmm::vect_norm2(rhs);
+ if (SOR_step) cout << "SOR step residual = " << precond_res << endl;
+ if (precond_res < residual) break;
+ cout << "Precond residual = " << precond_res << endl;
}
iter2.init();
// solving linear system for the global Newton method
if (0) {
- NS.compute_tangent_matrix(M, u);
- add_schwarz_mat<Matrixt, MatrixBi, Precond, local_solver>
- ASM(M, NS.get_vB(), iter4, P, iter.get_resmax());
- AS_global_solve(global_solver(), ASM, d, rhs, iter2);
+ NS.compute_tangent_matrix(M, u);
+ add_schwarz_mat<Matrixt, MatrixBi, Precond, local_solver>
+ ASM(M, NS.get_vB(), iter4, P, iter.get_resmax());
+ AS_global_solve(global_solver(), ASM, d, rhs, iter2);
}
else { // for exact gradient
- AS_exact_gradient<Matrixt, MatrixBi> eg(NS.get_vB());
- for (size_type i = 0; i < NS.get_vB().size(); ++i) {
- NS.compute_tangent_matrix(eg.vM[i], vx[i]);
- }
- eg.init();
- gmres(eg, d, rhs, gmm::identity_matrix(), 50, iter2);
+ AS_exact_gradient<Matrixt, MatrixBi> eg(NS.get_vB());
+ for (size_type i = 0; i < NS.get_vB().size(); ++i) {
+ NS.compute_tangent_matrix(eg.vM[i], vx[i]);
+ }
+ eg.init();
+ gmres(eg, d, rhs, gmm::identity_matrix(), 50, iter2);
}
// gmm::add(gmm::scaled(rhs, 0.1), u); ++iter;
external_ls.init_search(act_res, iter.get_iteration());
do {
- alpha = external_ls.next_try();
- gmm::add(gmm::scaled(d, alpha), u, x);
- NS.compute_F(rhs, x);
- act_res_new = gmm::vect_norm2(rhs);
+ alpha = external_ls.next_try();
+ gmm::add(gmm::scaled(d, alpha), u, x);
+ NS.compute_F(rhs, x);
+ act_res_new = gmm::vect_norm2(rhs);
} while (!external_ls.is_converged(act_res_new));
-
+
if (alpha != external_ls.converged_value()) {
- alpha = external_ls.converged_value();
- gmm::add(gmm::scaled(d, alpha), u, x);
- NS.compute_F(rhs, x);
- act_res_new = gmm::vect_norm2(rhs);
+ alpha = external_ls.converged_value();
+ gmm::add(gmm::scaled(d, alpha), u, x);
+ NS.compute_F(rhs, x);
+ act_res_new = gmm::vect_norm2(rhs);
}
if (iter.get_noisy() > 1) cout << endl;
- act_res = act_res_new;
- if (iter.get_noisy()) cout << "(step=" << alpha << ")\t unprecond res =
" << act_res << " ";
-
-
+ act_res = act_res_new;
+ if (iter.get_noisy())
+ cout << "(step=" << alpha
+ << ")\t unprecond res = " << act_res << " ";
+
++iter; gmm::copy(x, u);
}
}
