Télécharger la plaquette Cast3M


fr en

Télécharger arsluc.eso

Retour à la liste

Numérotation des lignes : 

arsluc

  1. C ARSLUC    SOURCE    CB215821  20/11/25    13:18:26     10792          
  2.       SUBROUTINE ARSLUC (IPRTRA,SIGMA,IPMAUP,QUAD,EPSI,INVER,
  3.      &                                  IPLVAR,IPLVER,IPLVAI,IPLVEI)
  4.  
  5. ************************************************************************
  6. *
  7. *                                A R S L U C
  8. *
  9. *
  10. *
  11. * FONCTION:
  12. * ---------
  13. *
  14. *     CONSTRUCTION DE LISTREEL'S SOLUTION DES FREQUENCES (REELLES ET/OU
  15. *     COMPLEXES) ET D'UN MLCHPO'S SOLUTION DES MODES CALCULES (REELS
  16. *     ET/OU COMPLEXES) A PARTIR DES VARIABLES DE SORTIE ARPACK
  17. *     (PROBLEME NON SYMETRIQUE)
  18. *
  19. *
  20. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE
  21. * -----------
  22. *
  23. *     IPRTRA  ENTIER     (E)    POINTEUR DES OPERATEURS DE TRAVAIL
  24. *
  25. *     SIGMA   COMPLEX DP (E)    VALEUR PROPRE DE SHIFT
  26. *
  27. *     INVER   LOGIQUE    (E)   .TRUE. -> PRODUIT SCALAIRE X'KX
  28. *                              .FALSE. -> PRODUIT SCALAIRE X'MX
  29. *
  30. *     IPMAUP  ENTIER     (E)    POINTEUR DES VARIABLES ARPACK
  31. *
  32. *     QUAD    LOGIQUE    (E)    PROBLEME QUADRATIQUE OU NON
  33. *
  34. *     IPLVAR  ENTIER     (S) POINTEUR DE L'OBJET 'LISTREEL' CONTENANT
  35. *                            LA SUITE DES FREQUENCES PROPRES RELLES
  36. *
  37. *     IPLVER  ENTIER     (S) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  38. *                            LA SUITE DES MODES PROPRES REELS
  39. *
  40. *     IPLVAI  ENTIER    (S)  POINTEUR DE L'OBJET 'LISTREEL' CONTENANT
  41. *                            LA SUITE DES FREQUENCES PROPRES IMAGINAIRES
  42. *
  43. *     IPLVEI  ENTIER    (S)  POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  44. *                            LA SUITE DES MODES PROPRES IMAGINAIRES
  45. *
  46. * SOUS-PROGRAMMES APPELES:
  47. * ------------------------
  48. *
  49. *   MOTS1,MAXIM1,VCH1,ARPSHI,NORMA1,ARPVER,MUCHPO
  50. *
  51. *
  52. * AUTEUR, DATE DE CREATION:
  53. * -------------------------
  54. *
  55. *     PASCAL BOUDA      11 JUILLET 2015
  56. *
  57. ************************************************************************
  58.  
  59.       IMPLICIT INTEGER(I-N)
  60.       IMPLICIT REAL*8 (A-H,O-Z)
  61.  
  62.  
  63. -INC PPARAM
  64. -INC CCOPTIO
  65. -INC CCREEL
  66. -INC SMRIGID
  67. -INC SMVECTD
  68. -INC SMLCHPO
  69. -INC SMLREEL
  70. -INC TARWORK
  71.  
  72.  
  73.  
  74.       INTEGER IPRTRA
  75.       COMPLEX*16 SIGMA
  76.       LOGICAL INVER
  77.       INTEGER IPMEUP
  78.       LOGICAL QUAD
  79.       INTEGER IPLVAR
  80.       INTEGER IPLVER
  81.       INTEGER IPLVAI
  82.       INTEGER IPLVEI
  83.  
  84.  
  85.       INTEGER IPRIGI
  86.       INTEGER JG, N1,INC
  87.       INTEGER IPCHO
  88.       INTEGER IPVEC
  89.       INTEGER IPLMOT
  90.       CHARACTER*(LOCOMP) MOTCLE
  91.       INTEGER TYPRO
  92.       INTEGER IPMONO,IPMON1,IPMON2
  93.       INTEGER IPMOD1,IPMOD2,IPMO22
  94.       REAL*8 MAXVAL, MAXVA1, MAXVA2
  95.       INTEGER N
  96.       COMPLEX*16 VSHIFT,FSHIFT
  97.       COMPLEX*16 CONJ
  98.       COMPLEX*16 ZERO
  99.  
  100.  
  101.       ZERO  =CMPLX(0.D0,0.D0)
  102.       VSHIFT=CMPLX(0.D0,0.D0)
  103.       FSHIFT=CMPLX(0.D0,0.D0)
  104.       i=1
  105.  
  106.       MRITRA=IPRTRA
  107.       SEGACT MRITRA
  108.  
  109.       IPRIGI=RIGI(1)
  110.  
  111.       SEGDES MRITRA
  112.  
  113.  
  114.       MRIGID=IPRIGI
  115.       SEGACT MRIGID
  116.       IPCHO=ICHOLE
  117.       SEGDES MRIGID
  118.  
  119.       MAUP=IPMAUP
  120.       SEGACT MAUP
  121.  
  122.       IF (QUAD) THEN
  123.         N=v(/1)/2
  124.         INCR=N
  125.       ELSE
  126.         N=v(/1)
  127.         INCR=0
  128.       ENDIF
  129.  
  130.       JG=dr(/1)
  131.       SEGINI MLREEL
  132.       IPLVAR=MLREEL
  133.       SEGDES MLREEL
  134.  
  135.       SEGINI MLREEL
  136.       IPLVAI=MLREEL
  137.       SEGDES MLREEL
  138.  
  139.       N1=dr(/1)
  140.       SEGINI MLCHPO
  141.       IPLVER=MLCHPO
  142.       SEGDES MLCHPO
  143.  
  144.       SEGINI MLCHPO
  145.       IPLVEI=MLCHPO
  146.       SEGDES MLCHPO
  147.  
  148. *Dans le cas de frequences complexes, il se peut que la nev-ieme
  149. *frequence soit complexe; dans ce cas, ce se sont pas nev mais nev+1
  150. *frequences qui sont calculees. Ainsi dans le cas ou la nev-ieme
  151. *frequence est reelle, on ajuste la taille des objets solutions
  152.       IF ( ABS(di(nev+1)) .EQ. 0.D0 ) THEN
  153.         JG=JG-1
  154.         N1=N1-1
  155.         MLREEL=IPLVAR
  156.         SEGADJ MLREEL
  157.         MLREEL=IPLVAI
  158.         SEGADJ MLREEL
  159.         MLCHPO=IPLVER
  160.         SEGADJ MLCHPO
  161.         MLCHPO=IPLVEI
  162.         SEGADJ MLCHPO
  163.       ENDIF
  164.  
  165.       DO 100 k=1,nev
  166.  
  167. *Lecture des valeurs propres: Si la valeur de shift est reelle ou
  168. *nulle, ARPACK fournit directement les valeurs propres. Dans le cas
  169. *contraire, il faut les recalculer a partir des quotients de Rayeigh
  170. *
  171. *En notant lambda=a+ib la valeur propre, V=X+iY le vecteur propre
  172. *associe, et B la matrice pour le produit scalaire, on a :
  173. *               - a=X*B*X+Y*B*Y
  174. *               - b=X*B*Y-Y*B*X
  175. * 18/08/2015 : Le shift complexe n'est pas implemente
  176.  
  177.  
  178.         IF (SIGMA .EQ. ZERO .OR. AIMAG(SIGMA) .EQ. 0.D0) THEN
  179. *Cas "simple": pas besoin de calculer les quotients
  180.           VSHIFT=CMPLX(dr(i),di(i))
  181.         ELSE
  182.           WRITE(IOIMP,*) 'Shift complexe non implemente !'
  183.           CALL ERREUR(5)
  184.         ENDIF
  185.  
  186.         IF (IIMPI .GT. 2) THEN
  187.           WRITE(IOIMP,*) ' '
  188.           WRITE(IOIMP,*) '--- Mode n°',k,' : on lit a l indice',i,' ---'
  189.           WRITE(IOIMP,*) 'Valeur propre',VSHIFT
  190.         ENDIF
  191.  
  192.  
  193.  
  194. *Conversion de la valeur propre de shift en frequence de shift
  195.         IF (IIMPI.GE.3) WRITE(IOIMP,*) '*** APPEL A ARPSHI ***'
  196.         CALL ARPSHI (FSHIFT,VSHIFT,QUAD,2)
  197. *Les valeurs propres sont classes dans les tableaux dr
  198. *(partie reelle) et di (partie imaginaire). Le tableau de sortie des
  199. *vecteurs propres d'ARPACK prend la convention suivante: si le mode est
  200. *reel, il occupe une colonne du tableau. S'il est complexe, il occupe
  201. *deux colonnes (une pour la partie reelle et une autre pour la partie
  202. *imaginaire.
  203.  
  204.  
  205.  
  206.  
  207.         IF ( ABS(di(i)) .EQ. 0.D0 ) THEN
  208.  
  209.           IF (IIMPI.GE.3) WRITE(IOIMP,*) 'ARSLUC : *** Cas reel ***'
  210.  
  211. *********************
  212. *******Cas reel******
  213. *********************
  214.  
  215. *La valeur propre est reelle, la frequence est reelle ou imaginaire pure
  216.           MLREEL=IPLVAR
  217.           SEGACT MLREEL*MOD
  218.           PROG(i)=REAL(FSHIFT)
  219.           IPLVAR=MLREEL
  220.           SEGDES MLREEL
  221.  
  222.           MLREEL=IPLVAI
  223.           SEGACT MLREEL*MOD
  224.           PROG(i)=AIMAG(FSHIFT)
  225.           IPLVAI=MLREEL
  226.           SEGDES MLREEL
  227.  
  228. *Recuperation du vecteur propre reel pur et transfomation en
  229. chpoint
  230.           INC=N
  231.           SEGINI MVECTD
  232.           DO j=1,N
  233.             VECTBB(j)=v(j+INCR,i)
  234. *Tracking : affichage du vecteur propre ARPACK
  235.             IF(IIMPI .GT. 300) WRITE(IOIMP,*) j,':',VECTBB(j)
  236.           ENDDO
  237.  
  238.           IPVEC=MVECTD
  239.           SEGDES MVECTD
  240.  
  241.           CALL VCH1 (IPCHO,IPVEC,IPMOD1,IPRIGI)
  242.           SEGSUP MVECTD
  243.  
  244. *La normalisation est impossible si le chpoint est nul
  245.           CALL MOTS1 (IPLMOT,MOTCLE)
  246.           CALL MAXIM1 (IPMOD1,IPLMOT,MOTCLE,0,MAXVAL)
  247.           IF (MAXVAL .LE. XPETIT) THEN
  248.             IPMONO=IPMOD1
  249.           ELSE
  250.             CALL NORMA1 (IPMOD1,IPLMOT,MOTCLE,IPMONO)
  251.           ENDIF
  252.  
  253. *Stockage dans le mlchpo
  254.           MLCHPO=IPLVER
  255.           SEGACT MLCHPO*MOD
  256.           ICHPOI(i)=IPMONO
  257.           IPLVER=MLCHPO
  258.           SEGDES MLCHPO
  259.  
  260. *Creeation d'un chpoint nul pour la partie imaginaire
  261.           INC=N
  262.           SEGINI MVECTD
  263.           DO j=1,N
  264.             VECTBB(j)=0.D0
  265.           ENDDO
  266.           IPVEC=MVECTD
  267.           SEGDES MVECTD
  268.           CALL VCH1 (IPCHO,IPVEC,IPMOD2,IPRIGI)
  269.           SEGSUP MVECTD
  270.  
  271. *pas de normalisation pour un chpoint nul
  272.           IPMONO=IPMOD2
  273.  
  274.           MLCHPO=IPLVEI
  275.           SEGACT MLCHPO*MOD
  276.           ICHPOI(i)=IPMONO
  277.           IPLVEI=MLCHPO
  278.           SEGDES MLCHPO
  279.  
  280. *Calcul de la norme et du residu du mode
  281.           TYPRO=iparam(7)
  282.           CALL ARPVER (IPRTRA,TYPRO,i,QUAD,.FALSE.,EPSI,INVER,
  283.      &                                            IPMOD1,IPMOD2,VSHIFT)
  284.  
  285.  
  286.         ELSE
  287.  
  288.           IF (IIMPI.GE.3) WRITE(IOIMP,*) 'ARSLUC : *** Cas complexe ***'
  289. *************************
  290. *******Cas complexe******
  291. *************************
  292.  
  293. *Les valeurs propres sont complexes conjuguees. On travaille par paires
  294. *Note: On "conjugue" la partie reelle dans le cas d'un probleme
  295. *quadratique aux valeurs propres
  296.           IF (.NOT. QUAD) THEN
  297.  
  298.             MLREEL=IPLVAR
  299.             SEGACT MLREEL*MOD
  300.             PROG(i)=REAL(FSHIFT)
  301.             PROG(i+1)=REAL(FSHIFT)
  302.             IPLVAR=MLREEL
  303.             SEGDES MLREEL
  304.  
  305.             MLREEL=IPLVAI
  306.             SEGACT MLREEL*MOD
  307.             PROG(i)=AIMAG(FSHIFT)
  308.             PROG(i+1)=-AIMAG(FSHIFT)
  309.             IPLVAI=MLREEL
  310.             SEGDES MLREEL
  311.  
  312.           ELSE
  313.  
  314.             MLREEL=IPLVAR
  315.             SEGACT MLREEL*MOD
  316.             PROG(i)=REAL(FSHIFT)
  317.             PROG(i+1)=-REAL(FSHIFT)
  318.             IPLVAR=MLREEL
  319.             SEGDES MLREEL
  320.  
  321.             MLREEL=IPLVAI
  322.             SEGACT MLREEL*MOD
  323.             PROG(i)=AIMAG(FSHIFT)
  324.             PROG(i+1)=AIMAG(FSHIFT)
  325.             IPLVAI=MLREEL
  326.             SEGDES MLREEL
  327.  
  328.           ENDIF
  329.  
  330. *Recuperation de la partie reelle du vecteur propre. On travaille par
  331. *paires (partie reele identique pour les deux)
  332.           INC=N
  333.           SEGINI MVECTD
  334.           DO j=1,N
  335.             VECTBB(j)=v(j+INCR,i)
  336.  
  337. *Tracking : affichage de la partie superieure du vecteur propre ARPACK
  338.             IF (IIMPI .GT. 300) THEN
  339.               WRITE(IOIMP,*) J,': RELAMX', v(j,i), 'RECALC',
  340.      &        REAL(VSHIFT)*v(j+INCR,i)-AIMAG(VSHIFT)*
  341.      &        v(j+INCR,i+1),'IMLAMX', v(j,i+1),'IMCALC',
  342.      &        AIMAG(VSHIFT)*v(j+INCR,i)+REAL(VSHIFT)*v(j+INCR,i+1)
  343.             ENDIF
  344.  
  345.           ENDDO
  346.           IPVEC=MVECTD
  347.           SEGDES MVECTD
  348.  
  349.           IF (IIMPI.GE.3) WRITE(IOIMP,*) '*** APPEL A VCH1 ***'
  350.           CALL VCH1 (IPCHO,IPVEC,IPMOD1,IPRIGI)
  351.           SEGSUP MVECTD
  352.  
  353. *Recuperation de la partie imaginaire du vecteur propre. On travaille
  354. *par paires (partie imaginaires opposees)
  355.           INC=N
  356.           SEGINI MVECTD
  357.           DO j=1,N
  358.             VECTBB(j)=v(j+INCR,i+1)
  359.  
  360. *Tracking : affichage de la partie inferieure du vecteur propre ARPACK
  361.             IF (IIMPI .GT. 500) THEN
  362.               WRITE(IOIMP,*) J,': REX',v(j+INCR,i),'IMX',v(j+INCR,i+1)
  363.             ENDIF
  364.  
  365.           ENDDO
  366.           IPVEC=MVECTD
  367.           SEGDES MVECTD
  368.  
  369.           CALL VCH1 (IPCHO,IPVEC,IPMOD2,IPRIGI)
  370.           SEGSUP MVECTD
  371.  
  372.           CALL MUCHPO (IPMOD2,-1.D0,IPMO22,1)
  373.  
  374.  
  375. *On recupere le max de la paire des chpoints
  376.           CALL MOTS1 (IPLMOT,MOTCLE)
  377.           CALL MAXIM1 (IPMOD1,IPLMOT,MOTCLE,0,MAXVA1)
  378.  
  379.           CALL MOTS1 (IPLMOT,MOTCLE)
  380.           CALL MAXIM1 (IPMOD2,IPLMOT,MOTCLE,0,MAXVA2)
  381.  
  382.           MAXVAL=MAX(MAXVA1,MAXVA2)
  383.  
  384. *Division de la partie reelle par le max
  385.           IF (MAXVA1 .EQ. 0.D0) THEN
  386.             IPMONO=IPMOD1
  387.           ELSE
  388.             CALL MUCHPO (IPMOD1,MAXVAL,IPMONO,-1)
  389.           ENDIF
  390.  
  391. *Stockage dans le mlchpo
  392.           MLCHPO=IPLVER
  393.           SEGACT MLCHPO*MOD
  394.           ICHPOI(i)=IPMONO
  395.           ICHPOI(i+1)=IPMONO
  396.           IPLVER=MLCHPO
  397.           SEGDES MLCHPO
  398.  
  399. *Division de la partie imaginaire par le max
  400.           IF (MAXVA2 .EQ. 0.D0) THEN
  401.             IPMON1=IPMONO
  402.             IPMON2=IPMONO
  403.           ELSE
  404.             CALL MUCHPO (IPMOD2,MAXVAL,IPMON1,-1)
  405. *Construction de la partie imaginaire du conjugue
  406.             CALL MUCHPO (IPMOD2,-MAXVAL,IPMON2,-1)
  407.           ENDIF
  408.  
  409. *Stockage dans de mlchpo
  410.           MLCHPO=IPLVEI
  411.           SEGACT MLCHPO*MOD
  412.           ICHPOI(i)=IPMON1
  413.           IPLVEI=MLCHPO
  414.           SEGDES MLCHPO
  415.  
  416. *Stockage dans de mlchpo
  417.           MLCHPO=IPLVEI
  418.           SEGACT MLCHPO*MOD
  419.           ICHPOI(i+1)=IPMON2
  420.           IPLVEI=MLCHPO
  421.           SEGDES MLCHPO
  422.  
  423. *Calcul de la norme et du residu des modes
  424.           TYPRO=iparam(7)
  425.  
  426.           CALL ARPVER (IPRTRA,TYPRO,i,QUAD,.FALSE.,EPSI,INVER,
  427.      &                                            IPMOD1,IPMOD2,VSHIFT)
  428.  
  429.           CONJ=CMPLX(REAL(VSHIFT),-AIMAG(VSHIFT))
  430.           CALL ARPVER (IPRTRA,TYPRO,(i+1),QUAD,.FALSE.,EPSI,INVER,
  431.      &                                            IPMOD1,IPMO22,CONJ)
  432.  
  433.  
  434.         ENDIF
  435.  
  436. *Critere de sortie de boucle
  437.         IF ( (i .EQ. JG) .OR. (i .EQ. (JG-1) .AND. di(i) .NE. 0) ) THEN
  438.           GOTO 101
  439.         ENDIF
  440.  
  441. *L'incrementation differe suivant le type de solution (reelle ou
  442. *complexe)
  443.         IF (di(i) .EQ. 0.D0) THEN
  444.           i=i+1
  445.         ELSE
  446.           i=i+2
  447.         ENDIF
  448.  
  449. 100   CONTINUE
  450. 101   CONTINUE
  451.  
  452. c       SEGDES MAUP
  453.       SEGDES MRITRA
  454.  
  455.       END
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales