191 lines
6.9 KiB
C++
191 lines
6.9 KiB
C++
|
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
|
||
|
|
/* ******** *** SparseLib++ */
|
||
|
|
/* ******* ** *** *** *** v. 1.5c */
|
||
|
|
/* ***** *** ******** ******** */
|
||
|
|
/* ***** *** ******** ******** R. Pozo */
|
||
|
|
/* ** ******* *** ** *** *** K. Remington */
|
||
|
|
/* ******** ******** A. Lumsdaine */
|
||
|
|
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
|
||
|
|
/* */
|
||
|
|
/* */
|
||
|
|
/* SparseLib++ : Sparse Matrix Library */
|
||
|
|
/* */
|
||
|
|
/* National Institute of Standards and Technology */
|
||
|
|
/* University of Notre Dame */
|
||
|
|
/* Authors: R. Pozo, K. Remington, A. Lumsdaine */
|
||
|
|
/* */
|
||
|
|
/* NOTICE */
|
||
|
|
/* */
|
||
|
|
/* Permission to use, copy, modify, and distribute this software and */
|
||
|
|
/* its documentation for any purpose and without fee is hereby granted */
|
||
|
|
/* provided that the above notice appear in all copies and supporting */
|
||
|
|
/* documentation. */
|
||
|
|
/* */
|
||
|
|
/* Neither the Institutions (National Institute of Standards and Technology, */
|
||
|
|
/* University of Notre Dame) nor the Authors make any representations about */
|
||
|
|
/* the suitability of this software for any purpose. This software is */
|
||
|
|
/* provided ``as is'' without expressed or implied warranty. */
|
||
|
|
/* */
|
||
|
|
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
|
||
|
|
|
||
|
|
#include <stdlib.h>
|
||
|
|
#include "diagpre_double_GR.h" // modif GR : découle d'une classe virtuelle
|
||
|
|
// pour avoir une interface standard
|
||
|
|
|
||
|
|
#include "Mat_pleine.h"
|
||
|
|
#include "MatBand.h"
|
||
|
|
#include "Mat_creuse_CompCol.h"
|
||
|
|
#include "ConstMath.h"
|
||
|
|
#include "MathUtil.h"
|
||
|
|
|
||
|
|
DiagPreconditioner_double::DiagPreconditioner_double(const CompRow_Mat_double &C)
|
||
|
|
: diag_(C.dim(0))
|
||
|
|
{
|
||
|
|
int i;
|
||
|
|
if ((i = CopyInvDiagonals(C.dim(0), &C.row_ptr(0), &C.col_ind(0),
|
||
|
|
&C.val(0), &diag(0))) != 0) {
|
||
|
|
cerr << "Diagonal preconditioner failure.";
|
||
|
|
cerr << " Zero detected in element " << i << endl;
|
||
|
|
Sortie(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
DiagPreconditioner_double::DiagPreconditioner_double(const CompCol_Mat_double &C)
|
||
|
|
: diag_ (C.dim (0))
|
||
|
|
{
|
||
|
|
int i;
|
||
|
|
if ((i = CopyInvDiagonals(C.dim(0), &C.col_ptr(0), &C.row_ind(0),
|
||
|
|
&C.val(0), &diag(0))) != 0) {
|
||
|
|
cerr << "Diagonal preconditioner failure.";
|
||
|
|
cerr << " Zero detected in element " << i << endl;
|
||
|
|
Sortie(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// cas d'une matrice abstraite
|
||
|
|
DiagPreconditioner_double::DiagPreconditioner_double(const Mat_abstraite &A) :
|
||
|
|
diag_ (A.Nb_ligne())
|
||
|
|
{ int i;
|
||
|
|
switch (A.Type_matrice())
|
||
|
|
{ case CREUSE_COMPRESSEE_LIGNE :
|
||
|
|
{const CompRow_Mat_double &C = (CompRow_Mat_double &) A;
|
||
|
|
// ici il y aura une modif a faire comme pour les matrices compressées colonne
|
||
|
|
// il faut passer la matrice dérivée et non la matrice mère
|
||
|
|
cout << "\n problème dans DiagPreconditioner_double(.. "
|
||
|
|
<< "voir le source .cpp, il y a du travail a faire";
|
||
|
|
Sortie(1);
|
||
|
|
if ((i = CopyInvDiagonals(C.dim(0), &C.row_ptr(0), &C.col_ind(0),
|
||
|
|
&C.val(0), &diag(0))) != 0)
|
||
|
|
{ cerr << "Diagonal preconditioner failure.";
|
||
|
|
cerr << " Zero detected in element " << i << endl;
|
||
|
|
Sortie(1);
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
case CREUSE_COMPRESSEE_COLONNE :
|
||
|
|
{const Mat_creuse_CompCol &C = (Mat_creuse_CompCol &) A;
|
||
|
|
if ((i = CopyInvDiagonals(C.dim(0), &C.col_ptr(0), &C.row_ind(0),
|
||
|
|
&C.val(0), &diag(0))) != 0)
|
||
|
|
{ cerr << "Diagonal preconditioner failure.";
|
||
|
|
cerr << " Zero detected in element " << i << endl;
|
||
|
|
Sortie(1);
|
||
|
|
}
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
default :
|
||
|
|
{if ((i = CopyInvDiagonals(A)) != 0)
|
||
|
|
{ cerr << "Diagonal preconditioner failure.";
|
||
|
|
cerr << " Zero detected in element " << i << endl;
|
||
|
|
Sortie(1);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
};
|
||
|
|
};
|
||
|
|
|
||
|
|
int
|
||
|
|
DiagPreconditioner_double::CopyInvDiagonals
|
||
|
|
(int n, const int *pntr, const int *indx
|
||
|
|
,const double *sa, double *diag)
|
||
|
|
{
|
||
|
|
int i, j;
|
||
|
|
|
||
|
|
for (i = 0; i < n; i++)
|
||
|
|
diag[i] = 0;
|
||
|
|
|
||
|
|
/* Find the diagonal elements */
|
||
|
|
for (i = 0; i < n; i++) {
|
||
|
|
for (j = pntr[i]; j < pntr[i+1]; j++) {
|
||
|
|
if (indx[j] == i) {
|
||
|
|
if (sa[j] == 0)
|
||
|
|
return i;
|
||
|
|
diag[i] = 1. / sa[j];
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if (diag[i] == 0)
|
||
|
|
return -i;
|
||
|
|
}
|
||
|
|
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
|
||
|
|
// surcharge de la fonction pour les matrices carrées et/ou bande
|
||
|
|
int DiagPreconditioner_double::CopyInvDiagonals(const Mat_abstraite& A)
|
||
|
|
{ switch (A.Type_matrice())
|
||
|
|
{ case CARREE :
|
||
|
|
break;
|
||
|
|
case CARREE_LAPACK :
|
||
|
|
break;
|
||
|
|
case BANDE_SYMETRIQUE :
|
||
|
|
case CARREE_SYMETRIQUE : break;
|
||
|
|
case RECTANGLE : break;
|
||
|
|
case BANDE_NON_SYMETRIQUE: break;
|
||
|
|
case CARREE_SYMETRIQUE_LAPACK : break;
|
||
|
|
case RECTANGLE_LAPACK : break;
|
||
|
|
case BANDE_SYMETRIQUE_LAPACK : break;
|
||
|
|
case BANDE_NON_SYMETRIQUE_LAPACK : break;
|
||
|
|
case CREUSE_NON_COMPRESSEE : break;
|
||
|
|
case CREUSE_COMPRESSEE_COLONNE : break;
|
||
|
|
case CREUSE_COMPRESSEE_LIGNE : break;
|
||
|
|
case DIAGONALE : break;
|
||
|
|
default :
|
||
|
|
cout << "\nErreur : valeur incorrecte du type de matrice !"
|
||
|
|
<< " = " << A.Type_matrice() << " , ce tye n'est pas encore "
|
||
|
|
<< " pris en compte dans le préconditionnement diagonal\n";
|
||
|
|
cout << "\n CopyInvDiagonals(Mat_abstraite& A, double *diag)\n";
|
||
|
|
Sortie(1);
|
||
|
|
};
|
||
|
|
int n = A.Nb_ligne();
|
||
|
|
for (int i = 1; i <= n; i++)
|
||
|
|
if (Abs(A(i,i)) < ConstMath::trespetit)
|
||
|
|
return -1;
|
||
|
|
else
|
||
|
|
diag(i-1) = 1./A(i,i);
|
||
|
|
return 0;
|
||
|
|
};
|
||
|
|
|
||
|
|
|
||
|
|
VECTOR_double
|
||
|
|
DiagPreconditioner_double::solve (const VECTOR_double &x) const
|
||
|
|
{
|
||
|
|
VECTOR_double y(x.size());
|
||
|
|
|
||
|
|
for (int i = 0; i < x.size(); i++)
|
||
|
|
y(i) = x(i) * diag(i);
|
||
|
|
|
||
|
|
return y;
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
VECTOR_double
|
||
|
|
DiagPreconditioner_double::trans_solve (const VECTOR_double &x) const
|
||
|
|
{
|
||
|
|
VECTOR_double y(x.size());
|
||
|
|
|
||
|
|
for (int i = 0; i < x.size(); i++)
|
||
|
|
y(i) = x(i) * diag(i);
|
||
|
|
|
||
|
|
return y;
|
||
|
|
}
|