[updated 9.Jan.2006]
Librairie valpro > Fichier ortogs.f |
SUBROUTINE ORTOGS (N,TYPE,KDEB,KFIN,MCX,MCQ,MCR,RST,CST)
Auteur : C. Hazard (Fevrier 1989), F.Mahe & D.Martin (janvier 1997)
Derniere modification : D.Martin (3 Mai 1999)
Orthonormalisation des colonnes d'une matrice avec ou sans construction de
la transformation ( algorithme "Twice Is Enough" )
Ce sous-programme orthonormalise les colonnes KDEB a KFIN de la matrice
de vecteurs (propres) colonnes X, les colonnes 1 a KDEB-1 etant supposees
deja orthonormales.
Il construit la matrice R (triangulaire superieure) de la transformation:
X = Q * R avec Q* Q = I
celle-ci n'est retournee a l'adresse MCR seulement si MCR.gt.0,
Il retourne la matrice Q orthogonale (qui peut coincider avec X si MCQ=MCX)
Il utilise la methode Gram-Schmidt dans un algorithme "twice is enough":
l'orthogonalisation est repetee lorsqu'une colonne de X voit sa norme reduite
d'un facteur superieur a COEF lors de la premiere orthogonalisation
(ici COEF=10).
Attention, l'algorithme ne permet pas necessairement de detecter une colonne
qui serait "presque" une combinaison lineaire des autres (voir la reference
pour les conditions assurees par l'algorithme).
Bibliographie : Parlett, The Symmetric Eigenvalue Problem,
Prentice-Hall, 1980, chap. 6, page 98
-- Arguments --
N nombre de lignes de la matrice des vecteurs colonnes X
TYPE type ('REEL' ou 'COMPLEXE') de la matrice
KDEB premiere colonne de la matrice a orthonormaliser
KFIN derniere colonne de la matrice a orthonormaliser
= nombre de colonnes de la matrice X
= nombre de lignes de la matrice carree triangulaire superieure R
MCX adresse de la matrice X
les colonnes du bloc (1:N,1:KFIN) de la matrice X donnees en entree
sont othonormalisees, et le resultat est range dans le bloc analogue
de Q en sortie.
Seules les colonnes KDEB a KFIN sont susceptibles d'etre modifiees,
les colonnes 1 a KDEB-1 etant supposees orthonormees en entree.
MCQ adresse de la matrice Q (coincide avec X si MCQ=MCX).
MCR adresse de la matrice R
en sortie le bloc (1:KFIN,KDEB:KFIN) de la matrice R contient la
transformation produite par l'orthonormalisation (matrice triangulaire
superieure).
Si lors de la seconde orthogonalisation, une colonne de X voit sa
norme reduite a nouveau d'un facteur superieur a COEF, l'element
diagonal de R correspondant et la colonne de X sont mis a 0.
CHARACTER*(*) TYPE
INTEGER N,KDEB,KFIN,MCX,MCQ,MCR
REAL RST(*)
COMPLEX CST(*)
INTEGER MCRDEB,MCRK,MCR1,MCQDEB,MCQK,MCQI,NORTHO,IST,K,I
REAL S,T,COEF2
COMPLEX C
COEF2=100.
IF (MCQ.NE.MCX)
& CALL TDANST (N*KFIN,TYPE,MCX,IST,RST,CST,TYPE,MCQ,IST,RST,CST)
MCRDEB=MCR+(KDEB-1)*KFIN
MCQDEB=MCQ+(KDEB-1)*N
adresse du 1er coeff. de la colonne KDEB de la matrice triang. sup. R
MCRK=MCRDEB
adresse du 1er coeff. de la colonne KDEB de la matrice orthogonale Q
MCQK=MCQDEB
DO 10 K=KDEB,KFIN
MCR1=MCRK
numero de l'orthogonalisation en cours
NORTHO=0
orthogonalisation de la colonne K / colonnes precedentes
1 NORTHO=NORTHO+1
T=0.
MCQI=MCQ
DO 3 I=1,K-1
produit scalaire des colonnes I et K
CALL TSCALT (N,TYPE,MCQI,RST,CST,TYPE,MCQK,RST,CST,C)
S=REAL(C)
T=T+S*S
mise a jour de la colonne K de X
CALL TPLUAT (N,TYPE,MCQK,RST,CST,'REEL',1,IST,-S,CST
& ,TYPE,MCQI,RST,CST,TYPE,MCQK,RST,CST)
MCQI=MCQI+N
remplissage de la matrice R
IF (MCR.GT.0.AND.NORTHO.EQ.1) THEN
RST(MCR1)=S
MCR1=MCR1+1
ENDIF
3 CONTINUE
carre de la norme de la nouvelle colonne K
CALL TSCALT (N,TYPE,MCQK,RST,CST,TYPE,MCQK,RST,CST,C)
S=REAL(C)
carre de la norme de l'ancienne colonne K
T=S+T
IF (S.LE.T/COEF2) THEN
IF (NORTHO.EQ.1) GOTO 1
IF (NORTHO.EQ.2) S=0.
ENDIF
normalisation
S=SQRT(S)
IF (MCR.GT.0) RST(MCRK+K-1)=S
IF (S.NE.0.) S=1./S
CALL TEQAT (N,TYPE,MCQK,RST,CST,'REEL',1,IST,S,CST
& ,TYPE,MCQK,RST,CST)
MCRK=MCRK+KFIN
MCQK=MCQK+N
10 CONTINUE
END
ortogs est appelé dans