init
This commit is contained in:
commit
c76662981c
30
TP1/cgs.m
Executable file
30
TP1/cgs.m
Executable file
|
@ -0,0 +1,30 @@
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
% ENSEEIHT - 1SN - Calcul scientifique
|
||||||
|
% TP1 - Orthogonalisation de Gram-Schmidt
|
||||||
|
% cgs.m
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
|
||||||
|
function Q = cgs(A)
|
||||||
|
|
||||||
|
% Recuperation du nombre de colonnes de A
|
||||||
|
[~, m] = size(A);
|
||||||
|
|
||||||
|
% Initialisation de la matrice Q avec la matrice A
|
||||||
|
Q = A;
|
||||||
|
|
||||||
|
% algo
|
||||||
|
|
||||||
|
Q(:,1) = Q(:,1) / norm(Q(:,1));
|
||||||
|
|
||||||
|
for j = 2:m
|
||||||
|
|
||||||
|
%projs = A(:,1:j-1)*A(:,1:j-1).'*A(:,j);
|
||||||
|
projs = Q(:,1:j-1)*Q(:,1:j-1).'*A(:,j);
|
||||||
|
|
||||||
|
Q(:,j) = Q(:,j) - projs;
|
||||||
|
|
||||||
|
Q(:,j) = Q(:,j) / norm(Q(:,j));
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
31
TP1/mgs.m
Executable file
31
TP1/mgs.m
Executable file
|
@ -0,0 +1,31 @@
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
% ENSEEIHT - 1SN - Calcul scientifique
|
||||||
|
% TP1 - Orthogonalisation de Gram-Schmidt
|
||||||
|
% mgs.m
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
|
||||||
|
function Q = mgs(A)
|
||||||
|
|
||||||
|
% Recuperation du nombre de colonnes de A
|
||||||
|
[~, m] = size(A);
|
||||||
|
|
||||||
|
% Initialisation de la matrice Q avec la matrice A
|
||||||
|
Q = A;
|
||||||
|
|
||||||
|
% Algo
|
||||||
|
|
||||||
|
for j = 1:m
|
||||||
|
|
||||||
|
for k = 1:j-1
|
||||||
|
|
||||||
|
proj = Q(:,k).'*Q(:,j)*Q(:,k);
|
||||||
|
|
||||||
|
Q(:,j) = Q(:,j) - proj;
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
Q(:,j) = Q(:,j) / norm(Q(:,j));
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
87
TP1/tp.m
Executable file
87
TP1/tp.m
Executable file
|
@ -0,0 +1,87 @@
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
% ENSEEIHT - 1SN - Calcul scientifique
|
||||||
|
% TP1 - Orthogonalisation de Gram-Schmidt
|
||||||
|
% tp.m
|
||||||
|
%--------------------------------------------------------------------------
|
||||||
|
|
||||||
|
clear
|
||||||
|
close all
|
||||||
|
clc
|
||||||
|
|
||||||
|
taille_ecran = get(0,'ScreenSize');
|
||||||
|
L = taille_ecran(3);
|
||||||
|
H = taille_ecran(4);
|
||||||
|
|
||||||
|
%% Calcul de la perte d'orthogonalite
|
||||||
|
|
||||||
|
% Rang de la matrice
|
||||||
|
n = 4;
|
||||||
|
|
||||||
|
% Puissance de 10 maximale pour le conditionnement
|
||||||
|
k = 16;
|
||||||
|
|
||||||
|
% Matrice U de test
|
||||||
|
U = gallery('orthog',n);
|
||||||
|
|
||||||
|
% Matrice de reference
|
||||||
|
D = eye(n);
|
||||||
|
|
||||||
|
% Initialisation de la matrice pour recuperer les pertes d'orthogonalite
|
||||||
|
po = zeros(2,k);
|
||||||
|
|
||||||
|
for i = 1:k
|
||||||
|
|
||||||
|
% Conditionnement de la matrice A
|
||||||
|
%TO DO: modifier D pour obtenir A tel K(A)=10^k
|
||||||
|
D(1,1) = 10^i;
|
||||||
|
A = U*D*U';
|
||||||
|
|
||||||
|
% Perte d'orthogonalite avec algorithme classique
|
||||||
|
Qcgs = cgs(A);
|
||||||
|
po(1,i) = norm(eye(n)-Qcgs'*Qcgs);
|
||||||
|
|
||||||
|
% Perte d'orthogonalite avec algorithme modifie
|
||||||
|
Qmgs = mgs(A);
|
||||||
|
po(2,i) = norm(eye(n)-Qmgs'*Qmgs);
|
||||||
|
|
||||||
|
|
||||||
|
% Perte d'orthogonalite avec algorithme classique x2
|
||||||
|
Qcgs2 = mgs(Qcgs);
|
||||||
|
po(3,i) = norm(eye(n)-Qcgs2'*Qcgs2);
|
||||||
|
|
||||||
|
% Perte d'orthogonalite avec algorithme modifié x2
|
||||||
|
Qmgs2 = mgs(Qmgs);
|
||||||
|
po(4,i) = norm(eye(n)-Qmgs2'*Qmgs2);
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
%% Affichage des courbes d'erreur
|
||||||
|
|
||||||
|
x = 10.^(1:k);
|
||||||
|
|
||||||
|
figure('Position',[0.1*L,0.1*H,0.8*L,0.75*H])
|
||||||
|
|
||||||
|
loglog(x,po(1,:),'r','lineWidth',2)
|
||||||
|
hold on
|
||||||
|
loglog(x,po(2,:),'b','lineWidth',2)
|
||||||
|
grid on
|
||||||
|
loglog(x,po(3,:),'g','lineWidth',2)
|
||||||
|
grid on
|
||||||
|
loglog(x,po(4,:),'y','lineWidth',2)
|
||||||
|
grid on
|
||||||
|
leg = legend('Gram-Schmidt classique',...
|
||||||
|
'Gram-Schmidt modifie',...
|
||||||
|
'Gram-Schmidt classique x2',...
|
||||||
|
'Gram-Schmidt modifie x2',...
|
||||||
|
'Location','NorthWest');
|
||||||
|
set(leg,'FontSize',14);
|
||||||
|
xlim([x(1) x(end)])
|
||||||
|
hx = xlabel('\textbf{Conditionnement $\mathbf{\kappa(A_k)}$}',...
|
||||||
|
'FontSize',14,'FontWeight','bold');
|
||||||
|
set(hx,'Interpreter','Latex')
|
||||||
|
hy = ylabel('$\mathbf{|| I - Q_k^\top Q_k ||}$','FontSize',14,'FontWeight','bold');
|
||||||
|
set(hy,'Interpreter','Latex')
|
||||||
|
title('Evolution de la perte d''orthogonalite en fonction du conditionnement',...
|
||||||
|
'FontSize',20)
|
||||||
|
|
||||||
|
|
29
TP2/Makefile
Normal file
29
TP2/Makefile
Normal file
|
@ -0,0 +1,29 @@
|
||||||
|
OBJLUPIV=TestLUPIV.o
|
||||||
|
LIBAUX = tlin/lib/libaux.a
|
||||||
|
F90 = gfortran
|
||||||
|
FOPTS = -g -fimplicit-none -ffixed-line-length-72 -fbounds-check -fbacktrace
|
||||||
|
|
||||||
|
all: TestLUPIV
|
||||||
|
|
||||||
|
TestLUPIV: laux $(OBJLUPIV)
|
||||||
|
@echo " .... Executable pour la facto LU AVEC pivotage : $@ "
|
||||||
|
$(F90) $(FOPTS) -o $@ $(OBJLUPIV) $(LIBAUX)
|
||||||
|
|
||||||
|
laux:
|
||||||
|
(cd tlin/lib; make)
|
||||||
|
|
||||||
|
clean:
|
||||||
|
(cd tlin/lib; make clean)
|
||||||
|
(cd tlin; rm -rf test *.o)
|
||||||
|
-rm TestLUPIV *.o
|
||||||
|
|
||||||
|
%.o: %.F90
|
||||||
|
$(F90) $(FOPTS) -c $<
|
||||||
|
|
||||||
|
%.o: %.f90
|
||||||
|
$(F90) $(FOPTS) -c $<
|
||||||
|
|
||||||
|
%.o: %.f
|
||||||
|
$(F90) $(FOPTS) -c $<
|
||||||
|
|
||||||
|
|
634
TP2/TestLUPIV.f90
Normal file
634
TP2/TestLUPIV.f90
Normal file
|
@ -0,0 +1,634 @@
|
||||||
|
! --------------------------------------------
|
||||||
|
! Programme principal permettant
|
||||||
|
! de tester la factorisation LU AVEC pivotage
|
||||||
|
! d'une matrice pleine
|
||||||
|
! --------------------------------------------
|
||||||
|
|
||||||
|
program mainlu
|
||||||
|
implicit none
|
||||||
|
interface
|
||||||
|
subroutine init(a, acopie, xexact, b, bcopie, x, y, r, &
|
||||||
|
p,q,lda, n,ipivot)
|
||||||
|
integer, intent(out) :: lda,n,ipivot
|
||||||
|
double precision, dimension(:,:), allocatable, intent(out) :: a, acopie
|
||||||
|
integer, dimension(:), allocatable , intent(out) :: p,q
|
||||||
|
double precision, dimension(:), allocatable, intent(out) :: &
|
||||||
|
xexact, b, bcopie, x, y, r
|
||||||
|
end subroutine init
|
||||||
|
end interface
|
||||||
|
|
||||||
|
! Dimension principale (taille de la structure de donnée allouee en memoire)
|
||||||
|
integer :: lda
|
||||||
|
! Dimension reelle du systeme a factoriser (on factorise A(1:n,1:n))
|
||||||
|
integer :: n
|
||||||
|
! Matrice a factoriser et permutation P et Q
|
||||||
|
double precision, dimension(:,:), allocatable :: a, acopie
|
||||||
|
integer, dimension (:), allocatable :: p, q
|
||||||
|
! Solution exacte, second membre, solution calculee du systeme et residu
|
||||||
|
double precision, dimension(:), allocatable :: xexact, b, x, &
|
||||||
|
r, bcopie
|
||||||
|
! Vecteur temporaire
|
||||||
|
double precision, dimension(:), allocatable :: y
|
||||||
|
! Determinant = DET_Mantisse * 2^DET_Exposant
|
||||||
|
! Pour eviter les underflows et overflows, le mantisse
|
||||||
|
! et l'exposant du determinant sont separees.
|
||||||
|
double precision :: det_mantisse
|
||||||
|
integer :: det_exposant
|
||||||
|
! Scalaire local
|
||||||
|
integer :: ierr
|
||||||
|
! Strat<EFBFBD>gie de pivot
|
||||||
|
integer :: ipivot
|
||||||
|
! EXTERNAL
|
||||||
|
external norme
|
||||||
|
double precision norme
|
||||||
|
!
|
||||||
|
! --------------------------------------------
|
||||||
|
! Initialisation matrice A et second membre b
|
||||||
|
! --------------------------------------------
|
||||||
|
! -- definir A et Xexact puis calculer b = A xexact
|
||||||
|
! (faire une copie de A --> Acopie et b --> bcopie)
|
||||||
|
! allouer x, y et r, P
|
||||||
|
|
||||||
|
call init(a, acopie, xexact, b, bcopie, x, y, r, p,q, lda, n,ipivot)
|
||||||
|
! -------------------------
|
||||||
|
! Factorisation de Gauss ( P A = L U ) et calcul du determinant
|
||||||
|
! -------------------------
|
||||||
|
!! write(6,*) ' '
|
||||||
|
!! write(6,*) ' ............................... '
|
||||||
|
!! write(6,*) ' FACTORISATION LU A ECRIRE '
|
||||||
|
!! write(6,*) ' ............................... '
|
||||||
|
call facto(a, lda, n,ipivot, p,q,ierr)
|
||||||
|
write(6,*) ' '
|
||||||
|
if(ierr==1)then
|
||||||
|
write(6,*) 'Echec de la factorisation'
|
||||||
|
|
||||||
|
else
|
||||||
|
! ----------------------
|
||||||
|
! Algorithme de descente ( L y = Pb)
|
||||||
|
! ----------------------
|
||||||
|
call descente ( a, lda, n, p, b, y)
|
||||||
|
! ----------------------
|
||||||
|
! Algorithme de remontee ( U x = y )
|
||||||
|
! ----------------------
|
||||||
|
call remontee ( a, lda, n, y, x)
|
||||||
|
|
||||||
|
!------------------------------------------
|
||||||
|
! Analyse de l'erreur :
|
||||||
|
! Calcul/affichage du residu || b - A x || et
|
||||||
|
! du residu equilibre || b - A x || / || |A||X| + |b| ||
|
||||||
|
! Calcul et affichage de || x - xexact<EFBFBD>|| / || xexact ||
|
||||||
|
!------------------------------------------
|
||||||
|
|
||||||
|
!! write(6,*) ' '
|
||||||
|
!! write(6,*) ' ............................... '
|
||||||
|
!! write(6,*) ' ROUTINE ANALYSE ERREUR A ECRIRE '
|
||||||
|
!! write(6,*) ' ............................... '
|
||||||
|
call analyse_erreur (acopie,lda,n,bcopie,x,r,xexact)
|
||||||
|
write(6,*) ' '
|
||||||
|
endif
|
||||||
|
|
||||||
|
deallocate (a, acopie, p, b, bcopie, x, r, y, xexact)
|
||||||
|
|
||||||
|
|
||||||
|
stop
|
||||||
|
end program mainlu
|
||||||
|
!--------------------------------------------------------------
|
||||||
|
!--------------------------------------------------------------
|
||||||
|
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
! Subroutines utilis<EFBFBD>es par le programme principal
|
||||||
|
!
|
||||||
|
!
|
||||||
|
!**************************************************
|
||||||
|
subroutine init(m,mcopie,xexact,b,bcopie, &
|
||||||
|
x,y,r,p,q,ldm,n,ipivot)
|
||||||
|
implicit none
|
||||||
|
! -------------------------------------------------------------
|
||||||
|
! Objectif:
|
||||||
|
! - Tous les tableaux/matrices passes en arguments
|
||||||
|
! sont alloues dans la routine.
|
||||||
|
! - Initialisation, de la matrice M, de la solution x
|
||||||
|
! et du vecteur b tel que b = M x
|
||||||
|
! M est une matrice de dimension principale ldm avec
|
||||||
|
! n colonnes.
|
||||||
|
! Deux classes de matrices sont proposees:
|
||||||
|
! Classe I/
|
||||||
|
! -Les matrices d'ordre pair sont diagonales dominantes
|
||||||
|
! -Les matrices d'ordre impair ont de petites valeurs
|
||||||
|
! sur la diagonale et peuvent conduire <EFBFBD> une perte
|
||||||
|
! de precision dans la solution
|
||||||
|
! Classe II/
|
||||||
|
! -12 matrices utilisant le generateur de LAPACK
|
||||||
|
! Soit x / x(i) = dble(i) on construit alors b = M x
|
||||||
|
! Une copie de M (-->Mcopie) et de b (--> bcopie) est effectuee
|
||||||
|
! --------------------------------------------------------------
|
||||||
|
!
|
||||||
|
interface
|
||||||
|
subroutine matgen(imat, a, n, info)
|
||||||
|
real(kind(1.d0)), allocatable :: a(:,:)
|
||||||
|
integer :: imat, n, info
|
||||||
|
end subroutine matgen
|
||||||
|
end interface
|
||||||
|
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(out) :: ldm,n,ipivot
|
||||||
|
double precision, allocatable, intent(out) :: &
|
||||||
|
m(:,:), mcopie(:,:), xexact(:), &
|
||||||
|
b(:), bcopie(:), x(:), y(:), r(:)
|
||||||
|
integer, allocatable, intent(out) :: p(:), q(:)
|
||||||
|
|
||||||
|
! Dimension maximale du systeme
|
||||||
|
integer :: ldmmax=5000
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
integer :: i,j, ierr, classmat, imat
|
||||||
|
logical :: pair
|
||||||
|
!
|
||||||
|
! Fonctions intrinseques
|
||||||
|
intrinsic dble, min
|
||||||
|
! -- taille du probleme
|
||||||
|
|
||||||
|
classmat = 1
|
||||||
|
ldm = -1
|
||||||
|
n = 0
|
||||||
|
do while (n.gt.ldm.or. ldm.ge.ldmmax)
|
||||||
|
write(6,*) ' Entrer la taille de la matrice lda (< ', &
|
||||||
|
ldmmax, ' ) ?'
|
||||||
|
read(5,*) ldm
|
||||||
|
write(6,*) ' Entrer la dimension reelle n (< ', ldm, ' ) ?'
|
||||||
|
read(5,*) n
|
||||||
|
!
|
||||||
|
write(6,*) ' Entrer la classe de matrice ', &
|
||||||
|
' (1 sinon 2 (generateur LAPACK)) ?'
|
||||||
|
read(5,*) classmat
|
||||||
|
!
|
||||||
|
if (classmat.eq.1) then
|
||||||
|
! -------------------------------
|
||||||
|
! Utilisation du generateur local
|
||||||
|
! -------------------------------
|
||||||
|
!
|
||||||
|
allocate(m(ldm,n), stat=ierr)
|
||||||
|
if (ierr > 0) then
|
||||||
|
write(6,*) ' Probleme allocation memoire :', &
|
||||||
|
' reduire n=', n
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
pair=(mod(n,2) .eq. 0)
|
||||||
|
! -- Initialisation de la matrice M
|
||||||
|
! Les matrices d'indice pair sont diagonales dominantes
|
||||||
|
! les matrices d'indice impair ont de petites valeurs
|
||||||
|
! sur la diagonale.
|
||||||
|
if (pair) then
|
||||||
|
do i = 1,n
|
||||||
|
do j=1,n
|
||||||
|
m(i,j) = dble(1)
|
||||||
|
end do
|
||||||
|
m(i,i) = dble(n+100)
|
||||||
|
end do
|
||||||
|
else
|
||||||
|
do i = 1,n
|
||||||
|
do j=1,n
|
||||||
|
m(i,j) = dble(i)*dble(j)
|
||||||
|
end do
|
||||||
|
m(i,i) = dble(i)*sqrt(epsilon(m(i,i)))
|
||||||
|
end do
|
||||||
|
endif
|
||||||
|
else
|
||||||
|
! -------------------------------------------------------------------
|
||||||
|
! Utilisation du generateur de matrices pour valider les codes LAPACK
|
||||||
|
! -------------------------------------------------------------------
|
||||||
|
imat = 0
|
||||||
|
do while (imat.le.0.or.imat.gt.12)
|
||||||
|
write(6,*) ' Generateur de matrices de LAPACK '
|
||||||
|
write(6,*) ' Entrer le numero de la matrice entre 1 et 12 ?'
|
||||||
|
read(5,*) imat
|
||||||
|
end do
|
||||||
|
! -- Initialisation de la matrice M
|
||||||
|
ldm = n
|
||||||
|
call matgen(imat, m, n, ierr)
|
||||||
|
if (ierr.ne.0) then
|
||||||
|
write(6,*) " Erreur dans le gnerateur LAPACK (matgen)"
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
! -- allocation des structures de données
|
||||||
|
allocate (p(n),q(n), mcopie(ldm,n), b(n), bcopie(n), x(n), &
|
||||||
|
r(n), y(n), xexact(n), stat=ierr)
|
||||||
|
|
||||||
|
! Stratégie de pivot
|
||||||
|
write(6,*) ' Entreer la strategie de pivot ', &
|
||||||
|
'0: sans pivot, 1 : pivot partiel, 2 : recherche du noyau'
|
||||||
|
read(5,*) ipivot
|
||||||
|
|
||||||
|
end do
|
||||||
|
if (ierr > 0) then
|
||||||
|
write(6,*) ' Probleme allocation memoire :', &
|
||||||
|
' reduire n=', n
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
! Copie de la matrice
|
||||||
|
mcopie = m
|
||||||
|
! -- Initialisation de la solution exacte dans xexact
|
||||||
|
do i=1,n
|
||||||
|
xexact(i) = dble(i)
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
! -- Calcul de b / b = A xexact
|
||||||
|
do i = 1, n
|
||||||
|
b(i) = dble(0)
|
||||||
|
do j=1,n
|
||||||
|
b(i) = b(i) + m(i,j) * xexact(j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
! Sauvegarde de b
|
||||||
|
bcopie = b
|
||||||
|
|
||||||
|
! -- impression des donnees
|
||||||
|
if (n.le.10) then
|
||||||
|
write (6,*) ' Matrice initiale A = '
|
||||||
|
do i = 1,min(n,6)
|
||||||
|
write(6,'(A,I3,6(E14.5))') ' ... Ligne ', i, m(i,1:min(n,6))
|
||||||
|
end do
|
||||||
|
write (6,'(A,6(E14.5))') ' xexact =', xexact(1:min(n,6))
|
||||||
|
write (6,'(A,6(E14.5))') ' b =', b(1:min(n,6))
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine init
|
||||||
|
!
|
||||||
|
!
|
||||||
|
!**************************************************
|
||||||
|
subroutine facto(m, ldm, n, ipivot, p, q, ierr)
|
||||||
|
implicit none
|
||||||
|
! ------------------------------------------------------
|
||||||
|
! Objectifs: Factorisation LU AVEC pivotage
|
||||||
|
! 1/ Calculer L / P M = L U
|
||||||
|
! La matrice M est de dimension principale ldm.
|
||||||
|
! On factorise le bloc M(1:n, 1:n)
|
||||||
|
! par l'algorithme de factorisation LU sans pivotage
|
||||||
|
! comme present<EFBFBD> en cours.
|
||||||
|
!
|
||||||
|
! Calculer le determinant
|
||||||
|
!
|
||||||
|
! En entr<EFBFBD>e la matrice M est <EFBFBD>cras<EFBFBD>e en sortie
|
||||||
|
! par les facteurs L et U.
|
||||||
|
!
|
||||||
|
! IERR > 0 si un pivot trop petit
|
||||||
|
! a <EFBFBD>t<EFBFBD> rencontr<EFBFBD> durant factorisation
|
||||||
|
! ------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(in) :: ldm,n,ipivot
|
||||||
|
double precision, intent(inout) :: m(ldm,n) ! matrice
|
||||||
|
integer, intent(out) :: p(n) ! permutation des lignes
|
||||||
|
integer, intent(out) :: q(n) ! permutation des colonnes
|
||||||
|
integer, intent(out) :: ierr
|
||||||
|
!
|
||||||
|
! Interface
|
||||||
|
interface
|
||||||
|
double precision function norme_f(a, n, m)
|
||||||
|
integer, intent(in) :: n,m
|
||||||
|
double precision, intent(in) :: a(n,m)
|
||||||
|
end function norme_f
|
||||||
|
end interface
|
||||||
|
|
||||||
|
! Fonctions intrinseques
|
||||||
|
intrinsic epsilon
|
||||||
|
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucle
|
||||||
|
integer :: i,j,k
|
||||||
|
integer :: t
|
||||||
|
|
||||||
|
double precision :: eps
|
||||||
|
|
||||||
|
write(6,*) " ========================="
|
||||||
|
write(6,*) " début de la FACTORISATION "
|
||||||
|
write(6,*) " ========================="
|
||||||
|
|
||||||
|
! epsilon machine (pour vérifier si pas d pivot trop petit)
|
||||||
|
eps = epsilon(eps)
|
||||||
|
|
||||||
|
! initialisation de la permutation (identité)
|
||||||
|
do k = 1, n
|
||||||
|
p(k) = k
|
||||||
|
end do
|
||||||
|
|
||||||
|
!! algo
|
||||||
|
do k = 1, n-1
|
||||||
|
|
||||||
|
select case (ipivot)
|
||||||
|
|
||||||
|
case (0)
|
||||||
|
if ( m(k, k) < eps * norme_f(m, n, n) ) stop ! on s'arrête si pivot trop petit
|
||||||
|
|
||||||
|
case (1)
|
||||||
|
i = maxloc(m(k:n,k), n-k) ! on cherche le pivot max
|
||||||
|
if ( m(i, k) < eps * norme_f(m, n, n) ) stop ! on s'arrête si pivot trop petit
|
||||||
|
|
||||||
|
! Swap rows i and k of A and b
|
||||||
|
t = p(k)
|
||||||
|
p(k) = p(i)
|
||||||
|
p(i) = t
|
||||||
|
|
||||||
|
case default
|
||||||
|
print*, "ipivot invalide"
|
||||||
|
stop
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
do i = k+1, n
|
||||||
|
m(i, k) = m(i, k) / m(k, k)
|
||||||
|
m(i, i) = m(i, i) - m(i, k) * m(k, i)
|
||||||
|
end do
|
||||||
|
|
||||||
|
end do
|
||||||
|
|
||||||
|
! -- impression des donnees :
|
||||||
|
! -- determinant
|
||||||
|
write(6,*) " =============================================="
|
||||||
|
write(6,*) " Fin de la FACTORISATION "
|
||||||
|
write(6,*) " ========================"
|
||||||
|
write(6,*) " ... AFFICHER DETERMINANT de la matrice "
|
||||||
|
|
||||||
|
! TODO afficher determinant
|
||||||
|
|
||||||
|
if (n.le.6) then
|
||||||
|
write (6,*) ' ... Matrice factorisee par ligne : '
|
||||||
|
do i = 1,n
|
||||||
|
write(6,'(6(E14.5))') m(i,1:n)
|
||||||
|
end do
|
||||||
|
write (6,*) " ... Permutation P : ", P(1:n)
|
||||||
|
write (6,*) " ... Permutation Q : ", Q(1:n)
|
||||||
|
endif
|
||||||
|
write(6,*) " =============================================="
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine facto
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
subroutine descente (m,ldm,n,p,b,x)
|
||||||
|
implicit none
|
||||||
|
! -----------------------------------------------------------
|
||||||
|
! Objectif:
|
||||||
|
! Resoudre L x = P b avec
|
||||||
|
! L la partie triangulaire inferieure de la matrice M,
|
||||||
|
! de dimension principale ldm et de n colonnes.
|
||||||
|
! Algorithme avec report:
|
||||||
|
! cet algorithme a assez naturellement tendance
|
||||||
|
! <EFBFBD> modifier le second membre ce qui n'est
|
||||||
|
! ni une bonne propriete ni
|
||||||
|
! indispensable comme le montre l'implementation fournie ici
|
||||||
|
! -------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(in) :: ldm,n
|
||||||
|
double precision, intent(in) :: m(ldm,n)
|
||||||
|
integer, intent(in) :: p(n)
|
||||||
|
! attention assez naturellement on va
|
||||||
|
! lors du report modifier b
|
||||||
|
! ( MAIS uniquement si on est peu attentif !)
|
||||||
|
double precision, intent(in) :: b(n)
|
||||||
|
double precision, intent(out) :: x(n)
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucle
|
||||||
|
integer :: i,j
|
||||||
|
double precision :: xtemp
|
||||||
|
!
|
||||||
|
! Algorithme avec report sans modification de b
|
||||||
|
!
|
||||||
|
|
||||||
|
x = b
|
||||||
|
|
||||||
|
! Ecrire la prise en compte de la permutation P en permutant
|
||||||
|
! le second membre qui est maintenant dans x
|
||||||
|
write(6,*)"Descente : permutation P non prise en compte"
|
||||||
|
|
||||||
|
|
||||||
|
do i = 1, n
|
||||||
|
! calcul du x(i)
|
||||||
|
! L a des 1 sur la diagonale (non stockes dans M!)
|
||||||
|
! report de la connaissance de x(i) sur es second membres.
|
||||||
|
! pour eviter de modifier b on note
|
||||||
|
! qu'a chaque etape i, les parties
|
||||||
|
! modifiees de b et calculees de x forment une partition de {1,..,n}
|
||||||
|
do j = i+1, n
|
||||||
|
! --- avec modification de b:
|
||||||
|
! --- b(j) = b(j) - M(j,i) * x(i)
|
||||||
|
! sans modification de b:
|
||||||
|
x(j) = x(j) - m(j,i) * x(i)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! write (6,*) " En sortie de la descente Ly = b: "
|
||||||
|
! write (6,'(A,6(E14.5))') ' b =', b(1:min(n,6))
|
||||||
|
! write (6,'(A,6(E14.5))') ' y =', x(1:min(n,6))
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine descente
|
||||||
|
!
|
||||||
|
!**************************************************
|
||||||
|
subroutine remontee (m,ldm,n,b,x)
|
||||||
|
implicit none
|
||||||
|
! ------------------------------------------------------
|
||||||
|
! Objectif:
|
||||||
|
! Resoudre U x = b avec
|
||||||
|
! U la partie triangulaire superieure de la matrice M,
|
||||||
|
! de dimension principale ldm et de n colonnes.
|
||||||
|
! Algorithme sans report.
|
||||||
|
! ------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(in) :: ldm,n
|
||||||
|
double precision, intent(in) :: m(ldm,n), b(n)
|
||||||
|
double precision, intent(out) :: x(n)
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucle
|
||||||
|
integer :: i,j
|
||||||
|
! Permutation
|
||||||
|
double precision :: xtemp
|
||||||
|
|
||||||
|
! Algorithme sans report.
|
||||||
|
do i = n, 1, -1
|
||||||
|
x(i) = b(i)
|
||||||
|
! on reporte la connaissance des x(j) pour j> i
|
||||||
|
! et on va chercher M(i,j) dans la triangulaire inferieure
|
||||||
|
do j = i+1, n
|
||||||
|
x(i) = x(i) - m(i,j) * x(j)
|
||||||
|
end do
|
||||||
|
x(i) = x(i) / m(i,i)
|
||||||
|
end do
|
||||||
|
|
||||||
|
! write (6,*) " En sortie de la remontee Ux = y: "
|
||||||
|
! write (6,'(A,6(E14.5))') ' y =', b(1:min(n,6))
|
||||||
|
! write (6,'(A,6(E14.5))') ' x =', x(1:min(n,6))
|
||||||
|
return
|
||||||
|
end subroutine remontee
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
subroutine analyse_erreur(m,ldm,n,b,x,r,xexact)
|
||||||
|
implicit none
|
||||||
|
! ------------------------------------------------------------
|
||||||
|
! Objectif:
|
||||||
|
! Analyse de l'erreur :
|
||||||
|
! Calcul et affichage de || b - M x || / || |M||X| + |b| ||
|
||||||
|
! Calcul et affichage de || x - xexact<EFBFBD>|| / || xexact ||
|
||||||
|
! M etant une matrice carree de
|
||||||
|
! dimension principale ldm et de n colonnes.
|
||||||
|
! -------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(in) :: ldm,n
|
||||||
|
double precision, intent(in) :: m(ldm,n), x(n), b(n)
|
||||||
|
double precision, intent(out) :: r(n)
|
||||||
|
double precision, intent(in) :: xexact(n)
|
||||||
|
|
||||||
|
! Interfaces
|
||||||
|
interface
|
||||||
|
double precision function norme ( x, n)
|
||||||
|
integer, intent(in) :: n
|
||||||
|
double precision, intent(in) :: x(n)
|
||||||
|
end function norme
|
||||||
|
|
||||||
|
double precision function norme_f ( a, n,m)
|
||||||
|
integer, intent(in) :: n,m
|
||||||
|
double precision, intent(in) :: a(n,m)
|
||||||
|
end function norme_f
|
||||||
|
end interface
|
||||||
|
|
||||||
|
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucles
|
||||||
|
integer :: i, j
|
||||||
|
double precision :: res_eq, err_rel
|
||||||
|
double precision, dimension(n) :: c, d
|
||||||
|
|
||||||
|
!!!!!!!!!!!!!!!!! A ECRIRE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
|
||||||
|
!! Calcul et affichage de des erreurs relatives directe et inverse
|
||||||
|
c = b - matmul(m, x)
|
||||||
|
d = matmul( abs(m), abs(x) ) + abs(b)
|
||||||
|
res_eq = sqrt( dot_product(c, c) / dot_product(d, d) )
|
||||||
|
!res_eq = norme_f(b - matmul(m,x), n, 1) / norme_f(matmul(abs(m),abs(x)) + abs(b), n, 1)
|
||||||
|
print*, "residu equilibré = ", res_eq
|
||||||
|
|
||||||
|
!! Calcul et afffichage de la norme du résidu
|
||||||
|
c = x - xexact
|
||||||
|
err_rel = sqrt( dot_product(c, c) / dot_product(xexact, xexact) )
|
||||||
|
print*, "erreur relative = ", err_rel
|
||||||
|
|
||||||
|
write(6,*) " "
|
||||||
|
write(6,*) " =============================================="
|
||||||
|
write(6,*) " Analyse d'ERREUR "
|
||||||
|
write(6,*) " ================ "
|
||||||
|
write(6,'(A,6(E14.5))') ' xexact =', xexact(1:min(n,6))
|
||||||
|
write(6,'(A,6(E14.5))') ' xcalcule =', x(1:min(n,6))
|
||||||
|
write(6,*) " =============================================="
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine analyse_erreur
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
subroutine noyau (m,ldm,n,q,kstop,x)
|
||||||
|
implicit none
|
||||||
|
! ------------------------------------------------------
|
||||||
|
! Objectif:
|
||||||
|
! Calculer une base du noyau depuis la factorisation incomplete
|
||||||
|
! Arret de la factosiation <EFBFBD> l'<EFBFBD>tape kstop
|
||||||
|
! Ker (AQ) = Im[-U11^(-1)*U12;Ir]
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! ------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Parametres :
|
||||||
|
! ----------
|
||||||
|
integer, intent(in) :: ldm,n,kstop
|
||||||
|
double precision, intent(in) :: m(ldm,n) ! facto LU incomplete de PAQ
|
||||||
|
integer, intent(in) :: q(n) ! permutation des colonnes
|
||||||
|
double precision, intent(out) :: x(n,n-kstop+1) ! base du noyau
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucle
|
||||||
|
integer :: i,j,l
|
||||||
|
|
||||||
|
!!!!!!!!!!!!!!!!! A ECRIRE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
!! Calcul d'une base du noyau de A
|
||||||
|
|
||||||
|
|
||||||
|
write (6,*) "Base du noyau non calculee"
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine noyau
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
double precision function norme (x, n)
|
||||||
|
implicit none
|
||||||
|
! ---------------------------------------
|
||||||
|
! Objectif: NORME = || x || , with 2 norm
|
||||||
|
! ----------------------------------------
|
||||||
|
integer, intent(in) :: n
|
||||||
|
double precision, intent(in) :: x(n)
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucles
|
||||||
|
integer :: i
|
||||||
|
!
|
||||||
|
! Fonctions intrinseques
|
||||||
|
intrinsic sqrt, dble
|
||||||
|
|
||||||
|
norme = dble(0)
|
||||||
|
do i=1, n
|
||||||
|
norme = norme + x(i) *x(i)
|
||||||
|
enddo
|
||||||
|
norme = sqrt(norme)
|
||||||
|
|
||||||
|
return
|
||||||
|
end function norme
|
||||||
|
|
||||||
|
!**************************************************
|
||||||
|
double precision function norme_f ( a, n,m)
|
||||||
|
implicit none
|
||||||
|
! ---------------------------------------
|
||||||
|
! Objectif: NORME_F = || A ||, avec norme de Frobenius
|
||||||
|
! ----------------------------------------
|
||||||
|
integer, intent(in) :: n, m
|
||||||
|
double precision, intent(in) :: a(n,m)
|
||||||
|
!
|
||||||
|
! Variables locales :
|
||||||
|
! -----------------
|
||||||
|
! Indices de boucles
|
||||||
|
integer :: i,j
|
||||||
|
!
|
||||||
|
! Fonctions intrinseques
|
||||||
|
intrinsic sqrt, dble
|
||||||
|
|
||||||
|
norme_f = dble(0)
|
||||||
|
do j=1,m
|
||||||
|
do i=1,n
|
||||||
|
norme_f = norme_f + a(i,j)*a(i,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
norme_f = sqrt(norme_f)
|
||||||
|
|
||||||
|
return
|
||||||
|
end function norme_f
|
||||||
|
|
||||||
|
|
||||||
|
|
Loading…
Reference in a new issue