[updated 3.Mar.2000]
*
SUBROUTINE ITINVA(NOMAT,NIVMAT,NOMM,NIVM,NOMMAC,NIVMAC,
* NOMVEC,NIVEC,SHIFT,ITMAX,REPSIL,VALPRO,NUIT)
*
* Auteur : F.Mahe 9 janvier 1992
* Derniere modification (F.Mahe) le 25 janvier 1996.
*
*
* NATURE DU CODE: Methode d'iteration inverse pour determiner la plus
* petite valeur propre et le vecteur propre normalise associe de
* la matrice de nom NOMAT et de niveau NIVMAT.
*
*
* DEMARCHE SUIVIE: Attention, l'acceleration de la convergence est ici une
* programmation dangeureuse car elle presuppose que certaines
* choses precises ont ete faites dans le programme principal
* et dans *Ptfixe*. Acceleration toutes les 7 iterations.
*
* PARAMETRES D'ENTREE
* NOMAT,NIVMAT nom et niveau de la matrice [A] (morse).
* sans donnee associee pour acceleration de la convergence(ASMTRM).
* ou construit selon les memes numerotations et type de stockage que
* NOMM,NIVM nom et niveau de la matrice de masse [M] (morse).
* NOMVEC,NIVEC nom et niveau du vecteur initial en entree
* remplace par le vecteur propre en sortie.
* NOMMAC,NIVMAC matrice 'ACSTTT' dep.du SHIFT pour acc.de la cvgce.
* SHIFT shift de selection de la valeur propre (modifie en sortie)
* ITMAX nombre d'iteration maxilal.
* REPSIL precision de l'iteration.
*
* PARAMETRES DE SORTIE
* VALPRO plus petite valeur propre (generalisee).
* NUIT nombre d'iterations.
*
*
CHARACTER*6 NOMAT,NOMM,NOMVEC,NOMMAC
CHARACTER TYPE,CARAC,CL*4,FORMEG3*10
INTEGER NIVMAT,NIVM,NIVEC,NUIT,IADRR,IAD,IBID
INTEGER NIVMAC
REAL REPSIL,VALPRO,KR,SHIFT
COMPLEX KC
*
INCLUDE 'CONTEX'
*
NIVESP=NIVESP+1
PREFIX(NIVESP)='*Itinva*'
*
* Algorithme d'iteration inverse *
*
* Vk = inv(A-SHIFT*I) Ukm1
* Uk = Vk/tvmax(Vk)
* a la convergence: tvmax(Vk) = + grande val.p.en module de inv(A-SHIFT*I)
* = 1/(+ petite val.p.en module de (A-SHIFT*I))
* = 1/((val.p.de A la + proche de SHIFT) - SHIFT)
*
*
* Normalisation initiale.
*
CALL MXTERM(NOMVEC,NIVEC,IADRR,KR,KC)
CALL PUTCST('f.m.norma','REEL',IBID,1/KR,KC,CARAC)
CALL AXTERM(NOMVEC,NIVEC,'f.m.norma','VVUKM1',1,-1)
*
* Decomposition 'LDtL' de [ADEUV]
*
CALL FALDLT(NOMAT,NIVMAT,'APROFI',1,-1) !syslin
*
* Recherche d'un vecteur propre, iteration.
*
NUIT=0 !compteur d'iterations.
100 NUIT=NUIT+1
*
* acceleration de la convergence : SHIFT <--- VALPRO.
* On suppose que dans le programme principal 'ACSTTT' de niveau 1 est definie
* par: 'ACSTTT' = 'RIGID' + '-K2'*'MASSE' + '-shift'*'SHIFTT'
* et que dans *Ptfixe* on a defini NOMAT:
* NOMAT = 'ADEUVV' = NOMMAC + [EFL], avec NOMMAC = 'ACSTTT'.
* Si pas de point fixe, NOMAT=NOMMAC='ADEUVV'.
*
IF (MOD(NUIT,7).EQ.0) THEN
* Calcul de la v.p. par produits scalaires:
* VALPRO = ([A]*(VECPRO), (VECPRO)) / ([M]*(VECPRO, (VECPRO))
CALL MATVEC(NOMAT,NIVMAT,'VVUKM1',1,'AVECPR',1,-1)
CALL MATVEC(NOMM,NIVM,'VVUKM1',1,'MVECPR',1,-1)
CALL XSTERM('S','AVECPR',1,'VVUKM1',1,TYPE,RPSA,KC)
CALL XSTERM('S','MVECPR',1,'VVUKM1',1,TYPE,RPSM,KC)
VALPRO=RPSA/RPSM
*
SHIFT=VALPRO+SHIFT
CALL PUTCST('-shift','REEL',IBID,-SHIFT,KC,CARAC)
* Pour appel dans Ptfixe: NOMMAC='ACSTTT' different de NOMAT.
IF ((NOMMAC(:6).NE.NOMAT(:6)).OR.(NIVMAC.NE.NIVMAT))
& CALL ASMTRM(NOMMAC,NIVMAC,-1)
CALL ASMTRM(NOMAT,NIVMAT,-1)
* Prise en compte des conditions essentielles eventuelles:
* Dirichlet homogene sur 'Gamma1', 'Gamma2', 'Gamma3' et 'Gamma4'.
CALL GETCST('CL',TYPE,IBID,KR,KC,CL)
IF (CL(1:1).EQ.'D')
& CALL CDESSE('PHIG1',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
IF (CL(2:2).EQ.'D')
& CALL CDESSE('PHIG2',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
IF (CL(3:3).EQ.'D') THEN
CALL GETCST('formeG3',TYPE,IBID,KR,KC,FORMEG3)
IF (FORMEG3(:8).EQ.'3parties') THEN
CALL CDESSE('PHIG3a',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
CALL CDESSE('PHIG3b',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
CALL CDESSE('PHIG3c',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
ELSE
CALL CDESSE('PHIG3',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
ENDIF
ENDIF
IF (CL(4:4).EQ.'D')
& CALL CDESSE('PHIG4',1,NOMAT,NIVMAT,NOMVEC,NIVEC)
*
CALL FALDLT(NOMAT,NIVMAT,'APROFI',1,-1) !syslin
ENDIF
*
* calcul du second membre: [M]Ukm1.
*
CALL MATVEC(NOMM,NIVM,'VVUKM1',1,'SECMEB',1,-1)
*
* Resolution de [A]Vk = [M]Uk-1
CALL FASV('APROFI',1,'SECMEB',1,NOMVEC,NIVEC,-1)
CALL GETERM(NOMVEC,NIVEC,IADRR,KR,KC)
VALPRO=1./KR
*
* Test d'arret
*
* nmax[Ukm1/tvmax(Vk)-VALPRO*Vk/tvmax(Vk)] < REPSIL
*
* car a la convergence Vk=Ukm1/VALPRO et tvmax(Vk) est un facteur
* de normalisation du test.
*
CALL MXTERM(NOMVEC,NIVEC,IADRR,KR,KC)
CALL PUTCST('f.m.norma','REEL',IBID,1/KR,KC,CARAC)
CALL PUTCST('facteur','REEL',IBID,-VALPRO/KR,KC,CARAC)
CALL CLTERM('VVUKM1',1,'f.m.norma',NOMVEC,NIVEC,
* 'facteur','VVUKM1',1,-1)
CALL MXTERM('VVUKM1',1,IBID,VALNI,KC)
VALNI=ABS(VALNI)
*
* Normalisation pour iteration suivante.
CALL AXTERM(NOMVEC,NIVEC,'f.m.norma','VVUKM1',1,-1)
IF (VALNI.GT.REPSIL) THEN
IF (NUIT.LT.ITMAX) THEN
GOTO 100
ENDIF
ENDIF
*
* Calcul de la v.p. par produits scalaires:
* VALPRO = ([A]*(VECPRO), (VECPRO)) / ([M]*(VECPRO, (VECPRO))
CALL MATVEC(NOMAT,NIVMAT,'VVUKM1',1,'AVECPR',1,-1)
CALL MATVEC(NOMM,NIVM,'VVUKM1',1,'MVECPR',1,-1)
CALL XSTERM('S','AVECPR',1,'VVUKM1',1,TYPE,RPSA,KC)
CALL XSTERM('S','MVECPR',1,'VVUKM1',1,TYPE,RPSM,KC)
VALPRO=RPSA/RPSM
*
* Test en norme 'max' de la convergence (verification):
* max([A]*(VECPRO) - VALPRO*[M]*(VECPRO)) < REPSIL
CALL PUTCST('facteur','REEL',IBID,-VALPRO,KC,CARAC)
CALL PUTCST('unite','REEL',IBID,1.,KC,CARAC)
CALL CLTERM('AVECPR',1,'unite','MVECPR',1,'facteur',
* NOMVEC,NIVEC,-1)
CALL MXTERM(NOMVEC,NIVEC,IAD,VALNI,KC)
VALNI=ABS(VALNI)
IF((VALNI.GT.REPSIL).AND.(NUIT.LT.ITMAX)) GOTO 100
*
* affectation des parametres de sorties
CALL T1TERM('VVUKM1',1,'=',NOMVEC,NIVEC,0)
VALPRO=VALPRO+SHIFT
*
* On ne tue pas les termes car ils seront utilises au appels
* suivants de *Itinv*: cela evite de generer des tableaux de
* compactage &plage et &PROFIL qui ne sont jamais detruits.
* CALL TBTUER('APROFI',1)
* CALL TBTUER('VVUKM1',1)
* CALL TBTUER('SECMEB',1)
* CALL TBTUER('AVECPR',1)
* CALL TBTUER('MVECPR',1)
*
NIVESP=NIVESP-1
*
RETURN
END
itinva est appelé dans (3 procédures)