C LANCZ1 SOURCE PV 11/03/08 21:15:43 6888 SUBROUTINE LANCZ1 (IPRIGI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) ************************************************************************ * * L A N C Z 1 * ----------- * * FONCTION: * --------- * * DETERMINATION D'UN VECTEUR DE LANCZOS. * * MODE D'APPEL: * ------------- * * CALL LANCZ1 (IPRIGI) * * PARAMETRES: (E)=ENTREE (S)=SORTIE * ----------- * * +A1 REEL DP (E) * +B1 REEL DP (E) * +IPMV1 ENTIER (E) * IPRIGI ENTIER (E) POINTEUR SUR LA MATRICE DE RIGIDITE * (OBJET 'RIGIDITE') EVENTUELLEMENT DECALEE. * +NUMV ENTIER (E) * +NUMV1 ENTIER (E) * +IPKVV ENTIER E/S LISTE COMPLETEE EN SORTIE. * CETTE LISTE EST SUPPRIMEE LORSQU'ELLE * DEVIENT INUTILE, C-A-D LORS DU DERNIER * APPEL A "LANCZ1". * +IPV1 ENTIER E/S NOUVEAU 'CHPOINT' EN SORTIE. * +IPV0 ENTIER (E) * (S) EGAL AU "IPV1" D'ENTREE. * +IPVTKV ENTIER E/S LISTE COMPLETEE EN SORTIE. * CETTE LISTE EST SUPPRIMEE LORSQU'ELLE * DEVIENT INUTILE, C-A-D LORS DU DERNIER * APPEL A "LANCZ1". * +IPVV ENTIER E/S LISTE COMPLETEE EN SORTIE. * +S0 REEL DP E/S * +S1 REEL DP E/S * * + = PARAMETRE PASSE PAR LE COMMUN "CLANCZ" (DEFINITIONS * DETAILLEES DANS LE SOUS-PROGRAMME "SIMUL0"). * * LEXIQUE: (ORDRE ALPHABETIQUE) * -------- * * COMBI2 SUBROUT. * IPV2 ENTIER POINTEUR SUR LE NOUVEAU VECTEUR DE LANCZOS, A * UN CERTAIN NIVEAU DE TRANSFORMATION. * IPV2A ENTIER ANALOGUE A "IPV2". * * MODE DE FONCTIONNEMENT: * ----------------------- * * CALCUL DU VECTEUR "V2" A PARTIR DE "V0" ET "V1" (QUI SONT LES 2 * VECTEURS DE LANCZOS PRECEDENTS): * * -1 * V2 = |K| .|M|.V1 - SIGNE(V1.|K|.V1).A.V1 - SIGNE(V0.|K|.V0).B.V0 * * AVEC: A = V1.|M|.V1 ET B = V0.|M|.V1 * |K| : MATRICE DE RIGIDITE, |M| : MATRICE MASSE. * * "V2" EST THEORIQUEMENT K-ORTHOGONAL A TOUS LES VECTEURS DE * LANCZOS CALCULES PRECEDEMMENT (DONT, NOTAMMENT, "V0" ET "V1"). * * SOUS-PROGRAMMES APPELES: * ------------------------ * * COMBI3, DTCHPO, DTLCHP, DTLREE, EXTRA4, INSER1, INSER4, LANCZ4, * MUCPRI, NORME2, ORTHO2 (OU ORTHO1 ?), RESOU1, YTMX. * * AUTEUR, DATE DE CREATION: * ------------------------- * * PASCAL MANIGOT 24 AVRIL 1985 * * LANGAGE: * -------- * * FORTRAN77 * ************************************************************************ * -INC CCOPTIO SEGMENT IDEMEN(0) * * REGOUPEMENT DE RENSEIGNEMENTS POUR LES SOUS-PROGRAMMES "LANCZ.": COMMON/CLANCZ/ A1,B1,S0,S1,NOMBRV,NUMV,NUMV1,IPV0,IPV1 & ,IPMV1,IPVV,IPKVV,IPVTKV,IPLMOX,IPLMOY * PARAMETER (UN = 1.D0 , XMOINS = -1.D0) PARAMETER (PRECIS = 1.D-15 , NFOIS = 1) * * * -- UNE ITERATION DE L'ALGORITHME DE LANCZOS -- * SEGINI IDEMEN IDEMEN(**)=IPMV1 NOID=0 NOEN=1 CALL RESOU1 (IPRIGI,IDEMEN,NOID,NOEN,1D-18) SEGACT IDEMEN IPV2 = IDEMEN(1) SEGSUP IDEMEN IF (IERR .NE. 0) RETURN CALL DTCHPO (IPMV1) * IF (NUMV .GT. 1) THEN * C IF (NUMV .EQ. 2) THEN C CALL LANCZ4 (IPV2,UN,IPV1,COEFF1, IPV2A) C ELSE C COEFF2 = XMOINS * B1 * S0 C CALL COMBI3 (IPV2,UN,IPV1,COEFF1,IPV0,COEFF2, IPV2A) C END IF C IF (IERR .NE. 0) RETURN COEFF1 = XMOINS * A1 * S1 CALL ADCHPO(IPV2,IPV1,IRETT,1.D0,COEFF1) IF (NUMV .EQ. 2) THEN IPV2A=IRETT ELSE COEFF2 = XMOINS * B1 * S0 CALL ADCHPO(IRETT,IPV0,IPV2A,1.D0,COEFF2) CALL DTCHPO(IRETT) END IF CALL DTCHPO (IPV2) IPV2 = IPV2A * IF (IIMPI .EQ. 369) THEN CALL NORME2 (IPV2,UN,IPRIGI, IPV2A,TEST1,IPLMOX, C IPLMOY,IPRX) IF (IERR .NE. 0) RETURN WRITE (IOIMP,2090) 2090 FORMAT (//' VERIFICATION, JUSTE APRES L''APPLICATION DE LA' & ,' FORMULE DE LANCZOS, DE L''ORTHOGONALITE DES '/ & ,' VECTEURS "VI" AVEC LE VECTEUR "V2" CREE (VALEURS ' & ,'A COMPARER A 1.) :'/) DO 80 IB80=1,NUMV1 CALL EXTRA4 (IPVV,IB80, IPVI) IF (IERR .NE. 0) RETURN CALL YTMX (IPVI,IPV2A,IPRIGI, TEST1) IF (IERR .NE. 0) RETURN WRITE (IOIMP,*) 'VI.K.V2 = ',TEST1 80 CONTINUE * END DO CALL DTCHPO (IPV2A) END IF * IPV0 = IPV1 S0 = S1 * END IF * * -- NORMALISATION ET, POUR COMPENSER LES ERREURS D'ARRONDI, * "K-ORTHOGONALISATION" -- * IF (NUMV .EQ. 1) THEN CALL NORME2 (IPV2,UN,IPRIGI, IPV1,S1,IPLMOX,IPLMOY,IPKV1) IF (IERR .NE. 0) RETURN CALL DTCHPO (IPV2) ELSE * HYP: NUMV > 1 CALL ORTHNG (IPV2,IPVV,IPVTKV,IPKVV,NUMV1,IPRIGI & ,S1,IPKV1, IPV1,IPLMOX,IPLMOY,precis,0) IF (IERR .NE. 0) RETURN CALL DTCHPO (IPV2) END IF * IF (IIMPI .EQ. 369) THEN WRITE (IOIMP,*) 'S1 = ',S1 END IF * ON REMPLACE "S1" (TRES PROCHE DE +/-1) PAR +/-1 S1 = SIGN(UN,S1) * CALL INSER4 (IPVV,NUMV,IPV1) IF (IERR .NE. 0) RETURN * CALL INSER1 (IPVTKV,NUMV,S1) IF (IERR .NE. 0) RETURN CALL INSER4 (IPKVV,NUMV,IPKV1) IF (IERR .NE. 0) RETURN END