C SESPA SOURCE BP208322 19/04/29 21:15:14 10213 C SESPA SOURCE WP 23//08/94 C SUBROUTINE SESPA ( IPLVAL, IPLVEC, NBMOD, IPRIGI, IPMASS ) ************************************************************************ * * SESPA * ----------- * * FONCTION: * --------- * * CONSTRUIT NBMOD ELEMENTS PROPRES EN ITERANT LE SOUS-ESPACE * IPLVEC, JUSQU'A LA CONVERGENCE DE CELUI-CI. * * MODE D'APPEL: * ------------- * * CALL SESPA ( IPLVAL, IPLVEC, NBMOD, IPRIGI, IPMASS ) * * PARAMETRES: (E)=ENTREE (S)=SORTIE * ----------- * * IPLVAL ENTIER (S) POINTEUR DE L'OBJET 'LISTREEL' CONTENANT * LA SUITE DE VALEURS PROPRES OBTENUES. * IPLVEC ENTIER (E) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT * LE SOUS-ESPACE INITIAL. * IPLVEC ENTIER (S) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT * LE SOUS-ESPACE FINAL. EN ORTHONORMALISANT * LES 'CHPOINT' DE CET ESPACE ON OBTIENT LES * VECTEURS PROPRES RECHERCHES. * NBMOD ENTIER (E) NOMBRE DE VECTEURS RECHERCHES. IPLVEC * CONTIENT PLUS QUE NBMOD 'CHPO', CAR CECI * PERMET DE CONVERGER PLUS RAPIDEMENT. * * IPRIGI ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' K * IPMASS ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' M * * * AUTEUR, DATE DE CREATION: * ------------------------- * * A.M. JOLIVALT, W. PASILLAS 13 / 07 / 94. ( ESOPE ) * *********************************************************************** SUBROUTINE SESPA ( IPLVAL, IPLVEC, NBMOD, IPRIGI, IPMASS ) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMLCHPO -INC SMLREEL ****** * -- CONSTANTES -- *** PARAMETER ( ITERMX = 40 ) ****** * -- ARGUMENTS -- *** POINTEUR IPLVEC.MLCHPO, IPLVAL.MLREEL INTEGER NBMOD, IPRIGI, IPMASS ****** * -- VARIABLES LOCALES -- *** POINTEUR IPLCH1.MLCHPO, IPLCH2.MLCHPO POINTEUR IPLVA1.MLREEL, IPLVA2.MLREEL INTEGER ILDIM LOGICAL BOOL **************************************************************** * INITIALISATIONS * **************************************************************** SEGACT ,IPLVEC ILDIM = IPLVEC.ICHPOI( /1 ) SEGDES ,IPLVEC ****** * -- ON PREND 0.0 COMME PREMIERE APPROXIMATION DES * VALEURS PROPRES -- *** JG = 0 SEGINI ,IPLVA1, IPLVA2 DO 90 IB90 = 1, ILDIM IPLVA1.PROG(**) = 0.D0 IPLVA2.PROG(**) = 0.D0 90 CONTINUE ****** * -- ON PREND LA LISTE IPLVEC COMME PREMIERE APPROX * DES MODES PROPRES -- *** N1 = 0 SEGINI ,IPLCH1 SEGACT ,IPLVEC DO 95 IB95 = 1, ILDIM IPCHPO = IPLVEC.ICHPOI(IB95) CALL COPIE2 ( IPCHPO, IPCH1 ) IF ( IERR .NE. 0 ) RETURN IPLCH1.ICHPOI(**) = IPCH1 95 CONTINUE SEGDES ,IPLVEC, IPLCH1 IPLCH2 = IPLVEC IPLVEC = 0 **************************************************************** * ITERATION DU SOUS-ESPACE IPLVEC * **************************************************************** ****** * -- REPETER JUSQU'A: * * CONVERGER. * * DEPASSER ITERMX ITERATIONS *** DO 100 IB100 = 1, ITERMX ****** * -- MISE A JOUR DE IPLVA1 -- *** CALL DTLREE ( IPLVA1 ) IF ( IERR .NE. 0 ) RETURN IPLVA1 = IPLVA2 ****** * -- MISE A JOUR DE IPLCH1 -- *** * -- ON DETRUIT -- SEGACT ,IPLCH1 DO 200 IB200 = 1, ILDIM IPCHPO = IPLCH1.ICHPOI(IB200) CALL DTCHPO (IPCHPO) IF ( IERR .NE. 0 ) RETURN 200 CONTINUE SEGDES ,IPLCH1 CALL DTLCHP ( IPLCH1 ) IF ( IERR .NE. 0 ) RETURN * -- ET ON RECOPIE -- SEGACT ,IPLCH2 N1 = 0 SEGINI ,IPLCH1 DO 300 IB300 = 1, ILDIM IPCHPO = IPLCH2.ICHPOI(IB300) CALL COPIE2 ( IPCHPO, IPCHP1 ) IF ( IERR .NE. 0 ) RETURN IPLCH1.ICHPOI(**) = IPCHP1 300 CONTINUE SEGDES ,IPLCH1, IPLCH2 ****** * -- UNE ITERATION DU SOUS-ESPACE IPLVEC -- *** CALL SESPA0 ( IPLCH2, IPLVA2, IPRIGI, IPMASS ) IF ( IERR .NE. 0 ) RETURN ****** * -- ON TESTE LA CONVERGENCE *** CALL SESPA5 ( IPLVA1, IPLVA2, IPLCH1, IPLCH2, 1 IPMASS, BOOL, NBMOD ) IF ( IERR .NE. 0 ) RETURN IF ( BOOL ) THEN * -- SI ON A CONVERGE, C'EST FINI ! -- IF ( IIMPI .EQ. 2 ) THEN WRITE ( IOIMP, 1000 ) IB100 1000 FORMAT( /1X, 'On a effectue ',I2,' iterations.', /) ENDIF GOTO 110 ENDIF IF ( IB100 .EQ. ITERMX ) THEN * -- SI NON, PAS DE CONVERGE, MAIS ON RENVOIE LA SOLUTION ! c WRITE ( IOIMP, 2000 ) ITERMX c 2000 FORMAT( /1X, 'Pas de convergence apres ',I2,' iterations.', c 1 /1X, 'La solution est quand meme renvoyee.', c 2 /1X, 'L''execution continue ...', / ) INTERR(1)=ITERMX CALL ERREUR(151) ENDIF 100 CONTINUE 110 CONTINUE ** estimation d'une borne superieure de l'erreur sur les valeurs propres ** (c.f. Argyris, Wilkinson) segact iplch2, iplva2, iplva1 do 50 ibmod = 1,nbmod ix = iplch2.ichpoi(ibmod) xlamda = IPLVA2.PROG( IBMOD) xlamd0 = IPLVA1.PROG( IBMOD) call mucpri(ix,ipmass,ixm1) call mucpri(ix,iprigi,ixk1) xco1= 1.d0 xco2 = -1d0*xlamda call adchpo(ixk1,ixm1,ires,xco1, xco2) call xtx1(ires,xres) call dtchpo(ixk1) call dtchpo(ixm1) call dtchpo(ires) dlamda = abs(xlamda - xlamd0) xerr1 = dlamda /abs(xlamda) xerr2 = (sqrt(abs(xres)))/abs(xlamda) IF ( (IIMPI .EQ. 2) .or. (.not. bool)) then reaerr(1) = xerr1 reaerr(2) = xerr2 interr(1) = ibmod c Valeur propre (omega**2) de rang ibmod : c convergence relative : xerr1 borne sup de l erreur relative : xerr2 call erreur(-327) write (ioimp,2010) 2010 format (1x) endif 50 continue segdes iplch2, iplva2,iplva1 **************************************************************** * NETTOYAGE DE LA MEMOIRE * **************************************************************** ****** * -- ON DETRUIT IPLVA1 ET IPLCH1 *** CALL DTLREE ( IPLVA1 ) IF ( IERR .NE. 0 ) RETURN SEGACT ,IPLCH1 DO 400 IB400 = 1, ILDIM IPCHPO = IPLCH1.ICHPOI(IB400) CALL DTCHPO (IPCHPO) IF ( IERR .NE. 0 ) RETURN 400 CONTINUE SEGDES ,IPLCH1 CALL DTLCHP ( IPLCH1 ) IF ( IERR .NE. 0 ) RETURN ****** * -- ON RENVOIE LES VALEURS ET VECTEURS PROPRES -- *** IPLVEC = IPLCH2 IPLVAL = IPLVA2 RETURN END