Télécharger la plaquette Cast3M


fr en

Télécharger vibrac.eso

Retour à la liste

Numérotation des lignes : 

vibrac

  1. C VIBRAC    SOURCE    BP208322  22/09/16    21:15:13     11454          
  2.       SUBROUTINE VIBRAC
  3. *
  4. ************************************************************************
  5. *      CALCUL DES MODES PROPRES DE PETITES MATRICES PLEINES            *
  6. *      MODES REELS OU COMPLEXES                                        *
  7. *      resp. PAR L'ALGORITHME DU QR ou QZ de Lapack                    *
  8. * _____________________________________________________________________*
  9. *                                                                      *
  10. *      AUTEUR :-VIBC calcul des modes Nicolas BENECH  le 24/02/95      *
  11. *              -option 'BALOU' Jerome CHARPENTIER le 29/11/96          *
  12. *              -amelioration rapidit� VIBC 'MODES' BP208322 16/12/2008 *
  13. *              -chercher des valeurs propres de la matrice monodromie  *
  14. *               LX236206   09/07/2014                                  *
  15. *              -suppression de l'option 'BALOU' BP, 09/09/2015         *
  16. *              -ajout nombre de modes souhaites BP, 13/01/2016         *
  17. *              -ajout calcul des vp reelles d'1 rigidite, BP, 2022-09  *
  18. * _____________________________________________________________________*
  19. *                                                                      *
  20. * APERCU DES SYNTAXES (details --> voir notice) :                      *
  21. *                                                                      *
  22. * SYNTAXE 1 :                                                          *
  23. *     BAS2 = VIBC  MASS1 RIG1 (AMOR1) (BAS1) (ENT1);                   *
  24. *                                                                      *
  25. *     RESOLUTION DE [K + (i*2*pi*w)*C - (2*pi*w)**2 M] X = 0           *
  26. *                                                                      *
  27. *     EN METTANT LE PB SOUS LA FORME D'ETAT :                          *
  28. *            [ A -  \lambda B ] . Y = 0                                *
  29. *                                                                      *
  30. *     OU     A = [ -K   0  ]      et    B = [ C   M ]                  *
  31. *                [  0   M  ]                [ M   0 ]                  *
  32. *                                                                      *
  33. * _____________________________________________________________________*
  34. *                                                                      *
  35. *  SYNTAXE 2 :                                                         *
  36. *      BAS2 = VIBC  RIG1 (RIG2 RIG3 RIG4) ;                            *
  37. *                                                                      *
  38. *      RESOLUTION DE [ A - \lambda I ] . X = 0                         *
  39. *                                                                      *
  40. *      OU   A = |   RIG1            si seul RIG1 est fourni            *
  41. *               |                                                      *
  42. *               | [ RIG1  RIG2 ]    si 4 rigidites sont fournies       *
  43. *               | [ RIG3  RIG4 ]                                       *
  44. *                                                                      *
  45. ************************************************************************
  46. *
  47. ************************************************************************
  48. *     DECLARATIONS
  49. ************************************************************************
  50.  
  51.       IMPLICIT INTEGER(I-N)
  52.       IMPLICIT REAL*8 (A-H,O-Z)
  53.  
  54.       DATA ZERO / 0.0D0 /      
  55.       INTEGER ERR, NUMAFF
  56.       LOGICAL getVEC, doRECO, AFFICH
  57.       REAL*8 XPREC1
  58.  
  59. -INC PPARAM
  60. -INC CCOPTIO
  61. -INC SMRIGID
  62.       POINTEUR MAT1.MRIGID, MAT2.MRIGID
  63.       POINTEUR CHPO1.MCHPOI, CHPO2.MCHPOI, CHPO4.MCHPOI
  64.       POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  65. *     tableaux pour Lapack
  66.       SEGMENT MWORK
  67.         REAL*8 LAMBDA(NWORK),WORK(LWORK)
  68.       ENDSEGMENT
  69.  
  70.  
  71. ************************************************************************
  72. *     INITIALISATIONS
  73. ************************************************************************
  74.  
  75. *     indique si l'on calcule ou non les vecteurs propres
  76.       getVEC = .TRUE.
  77. *     indique si l'on recombine ou non les vecteurs propres
  78.       doRECO = .FALSE.
  79. *     indique si l'on affiche ou non des messages de deroulement
  80. c       AFFICH = .FALSE.
  81.       AFFICH = (IIMPI.GE.1)
  82. *     indique la syntaxe 
  83. *     =0 : erreur
  84. *     =1 : syntaxe 1      [K + (i*2*pi*w)*C - (2*pi*w)**2 M] X = 0
  85. *     =21 : syntaxe 2.1   [K - lambda I] X = 0
  86. *     =22 : syntaxe 2.2   [[K1 K2 ; K3 K4] - lambda [I 0 ; 0 I]] X = 0 (matrice de monodromie par ex.)
  87.       ISYNT=0
  88.       IMONO=0
  89.       MWORK=0
  90.  
  91.  
  92. ************************************************************************
  93. *     LECTURE DES ARGUMENTS
  94. ************************************************************************
  95.  
  96. *     RIGIDITES : la syntaxe dépend du nombre de rigidités lues
  97.       CALL LIROBJ('RIGIDITE',IPMAS,1,IRET)
  98.       IF (IERR.NE.0) RETURN
  99.       CALL LIROBJ('RIGIDITE',IPRIG,0,IRET)
  100.       IF (IERR.NE.0) RETURN
  101. *   - n=1 matrice lue
  102.       IF(IRET.eq.0) THEN
  103.          IF (AFFICH) write(ioimp,*) 'syntaxe : [A - lambda*I] X = 0'
  104.          ISYNT=21
  105.          GOTO 1000
  106. *   - n>1 matrices lues
  107.       ELSE
  108. *   - n=2 ou 3 matrices lues (syntaxe par defaut)
  109.         ISYNT=1
  110.         CALL LIROBJ('RIGIDITE',IPAMO,0,IRET)
  111.         IF (IERR.NE.0) RETURN
  112.         IF (IRET.ne.0) then
  113.           CALL LIROBJ('RIGIDITE',IPMON,0,IMONO)
  114.           IF (IERR.NE.0) RETURN
  115. *   - n=4 matrices lues
  116.           IF (IMONO.eq.1) THEN
  117.             IF (AFFICH) write(ioimp,*) 'syntaxe : monodromie'
  118.             ISYNT=22
  119.             GOTO 1000
  120.           ENDIF
  121.         ENDIF
  122.       ENDIF    
  123.  
  124. *     ici, seulement synatxe 1 :
  125.       IF (AFFICH) write(ioimp,*) 'syntaxe : [K +iwC -w^2M] X = 0'
  126.  
  127. *     TABLE DE BASE MODALE
  128.       CALL LIROBJ('TABLE   ',IBASR,0,IRET)
  129.       IF (IERR.NE.0) RETURN
  130.       doRECO = (IRET.EQ.1)
  131.  
  132. *     NOMBRE DE MODES SOUHAITES
  133.       NWANTED=0
  134.       CALL LIRENT(NWANTED,0,IRET)
  135.       IF (IERR.NE.0) RETURN
  136.  
  137. *     TRI DES 3 (OU 2) RIGIDITES 
  138.       CALL QZTRIR(IPMAS, IPRIG, IPAMO)
  139.       GOTO 1000
  140.  
  141.  
  142.  1000 CONTINUE
  143. ************************************************************************
  144. *     ASEMBLAGE DES MATRICES A (ET B) 
  145. *     en fonction des syntaxes
  146. ************************************************************************
  147.       if (ISYNT.EQ.1) then   
  148.         CALL QZCONS(IPMAS,IPRIG,IPAMO,MATA,MATB,IPT1,MCOMP,.false.)
  149.         GOTO 3000
  150.       elseif (ISYNT.EQ.22) then
  151.         CALL QZCON3(IPMAS,IPRIG,IPAMO,IPMON,MATA,MATB,IPT1,MCOMP,.true.)
  152.         GOTO 3000
  153.       elseif (ISYNT.EQ.21) then
  154.         CALL QZCON1(IPMAS,MATA,MATB,IPT1,MCOMP)
  155.         GOTO 2000
  156.       else
  157.         CALL ERREUR(5)
  158.       endif
  159.       IF (IERR.NE.0) RETURN
  160.  
  161.  
  162.  2000 CONTINUE
  163. ************************************************************************
  164. *     ALGORITHME DU QR de LAPACK 
  165. *    (pourrait etre dans une subroutine nommée vibc1)
  166. ************************************************************************
  167.  
  168. *     CREATION SEGMENT DE TRAVAIL
  169.       NDDL=MATA.RE(/1)
  170. *     on dimensionne nb identiquement a DSYTRD
  171.       NB=ILAENV(1,'DSYTRD','U', NDDL, -1, -1, -1 )
  172.       NWORK=NDDL
  173.       LWORK=(NB+2)*NDDL
  174.       SEGINI,MWORK      
  175.  
  176. *     APPEL A DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
  177.       CALL DSYEV('V','U',NDDL,MATA.RE(1,1,1),NDDL,
  178.      &           LAMBDA,WORK,LWORK,IRET)
  179.  
  180. *     gestion erreur Lapack -> erreur Cast3M
  181.       IF (IRET.GT.0) THEN
  182.         CALL ERREUR(26)
  183.         RETURN
  184.       ELSEIF(IRET.LT.0) THEN
  185.         MOTERR(1:8)='VIBC    '
  186.         CALL ERREUR(997)
  187.         RETURN
  188.       ENDIF
  189.  
  190. *      write(*,*) 'eigenvalues =',(LAMBDA(iou),iou=1,NDDL)
  191. *      write(*,*) 'Vector #',(j,j=1,NDDL) 
  192. *      do i=1,NDDL
  193. *      write(*,*) (MATA.RE(i,j,1),j=1,NDDL)
  194. *      enddo
  195.  
  196. *     CREATION D'UNE BASE DE MODES PROPRES COMPLEXES POUR VIBC
  197.       CALL QRBASR (MWORK,MATA,IPT1,MCOMP,IBASC,NDDL)
  198.       IF (IERR.NE.0) RETURN
  199.  
  200. *     MENAGE
  201.       SEGSUP,MATA,MATB,MWORK
  202.  
  203. *     FIN NORMALE      
  204.       GOTO 9000
  205.  
  206.  
  207.  3000 CONTINUE
  208. ************************************************************************
  209. *     ALGORITHME DU QZ 
  210. *    (pourrait etre dans une subroutine nommée vibc2)
  211. ************************************************************************
  212.  
  213. * -------------------------------------------------------------------
  214. *     MISE SOUS FORME DE HESSENBERG SUPERIEURE DE MATA
  215. *     ET TRIANGULARISATION SUPERIEURE SIMULATANEE DE MATB
  216. *premier pas
  217.       CALL QZHESS(MATA,MATB,getVEC,IPRI3)
  218.       IF (IERR.NE.0) RETURN
  219.       XPREC1 = ZERO
  220.       CALL LIRREE(XPREC1,0,IRETOU)
  221. *deuxieme pas
  222.       CALL QZREDU(MATA,MATB,XPREC1,getVEC,IPRI3,ERR)
  223.       IF (IERR.NE.0) RETURN
  224.       IF (ERR.NE.0)  WRITE (*,*) 'Mode ',ERR,' : trop d iterations !'
  225.      &,' L execution continu ...'
  226.  
  227.  
  228.  
  229. * -------------------------------------------------------------------
  230. *     EXTRACTION DES VALEURS PROPRES
  231. *troisieme pas
  232.       CALL QZVALP(MATA,MATB,IPAR,IPAI,IPBE,getVEC,IPRI3)
  233.       IF (IERR.NE.0) RETURN
  234.  
  235. * -------------------------------------------------------------------
  236. *     CALCUL DES VECTEURS PROPRES COMPLEXES
  237. *quatrieme pas
  238.       CALL QZVECP(MATA,MATB,IPAR,IPAI,IPBE,IPRI3)
  239.       IF (IERR.NE.0) RETURN
  240.  
  241. * -------------------------------------------------------------------
  242. *     CREATION D'UNE BASE DE MODES PROPRES COMPLEXES POUR VIBC
  243.       CALL QZBASC(IPAR,IPAI,IPBE,IPRI3,IPT1,MCOMP,IBASC,ERR,NWANTED)
  244.       IF (IERR.NE.0) RETURN
  245.  
  246. * -------------------------------------------------------------------
  247. *     RECOMBINAISON D'UNE BASE DE MODES COMPLEXES A PARTIR
  248. *     D'UNE BASE DE MODES REELS (= RECOMBINAISON MODALE)
  249.       IF (doRECO) THEN
  250.         CALL QZREST(IBASR, IBASC)
  251.       ENDIF
  252.  
  253.  
  254.  
  255.  9000 CONTINUE
  256. ************************************************************************
  257. *     ECRITURE DU RESULTAT
  258. ************************************************************************
  259.  
  260.       IF (AFFICH) write(ioimp,*) 'ecriture de la table #',IBASC
  261.       CALL ECROBJ('TABLE   ',IBASC)
  262.       RETURN
  263.  
  264.  
  265.       RETURN
  266.       END
  267.  
  268.  
  269.  
  270.  

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