Sources Esope | Cast3M

exdiag
C EXDIAG    SOURCE    GOUNAND   25/04/30    21:15:04     12258          
      SUBROUTINE EXDIAG(IOPT)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C*************************************************************************
C     Operateur EXDIAG
C
C     OBJET   : Extrait une diagonale d'une matrice au format
C               MATRIK
C               Pour cela, on doit l'assembler.
C               L'assemblage se fait comme dans KRES2.
C               Si IOPT=1, on renvoie la diagonale de la matrice.
C               Si IOPT=2, on renvoie un inverse approché (SPAI)
C                 diagonal de la matrice entrée.
C                 Son expression est :
C                 m_ii = a_ii / \sum_{j=1,n} aij^2
C
C     Voir aussi la notice de KOPS
C
C***********************************************************************
C HISTORIQUE : 21/03/08 : Première version
C
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMCHPOI
      POINTEUR MCHDIA.MCHPOI
-INC SMMATRIK
      POINTEUR MATRI2.MATRIK
      POINTEUR AMORS.PMORS
      POINTEUR AISA.IZA
      POINTEUR IDIAGO.IZA
      POINTEUR ISCAR.IZA
      INTEGER IMPR,IRET
C
      CHARACTER*4 MRENU,MMULAG
      CHARACTER*8 TYPE
      CHARACTER*8 TYMATK,TYRIGI,BLAN
      DATA TYMATK,TYRIGI,BLAN/'MATRIK  ','RIGIDITE','        '/
*
      IMPR=0
*
* Lecture de la matrice
*
      TYPE=BLAN
      CALL QUETYP(TYPE,1,IRET)
      IF (IRET.EQ.0) GOTO 9999
      IF (TYPE.EQ.TYRIGI) THEN
         IF (IOPT.EQ.2) THEN
* Transformation en matrik en changement de noms d'inconnues
            CALL RIMA
            IF (IERR.NE.0) RETURN
            CALL MACHI2(1)
            IF (IERR.NE.0) RETURN
         ELSE
            CALL LIROBJ('RIGIDITE',IRIG,1,IRET)
            IF (IERR.NE.0) RETURN
*
            CALL EXDIAR(IRIG,ICHP)
            IF (IERR.NE.0) RETURN
*
            CALL ECROBJ('CHPOINT ',ICHP)
            RETURN
         ENDIF
      ENDIF
      TYPE=TYMATK
      CALL LIROBJ(TYPE,MATRIK,1,IRET)
      IF(IRET.EQ.0) GOTO 9999
C
C Assemblage proprement dit
C
C  Attention, on effectue une recopie du chapeau pour ne garder
C aucune information de préconditionnement (assemblage, numérotation)
C dans la matrice originale sinon une résolution subséquente poserait
C problème !!!!!!
C
      SEGINI,MATRI2=MATRIK
      MATRIK=MATRI2
      MATASS=MATRIK
      MRENU='RIEN'
      MMULAG='RIEN'
*     SG 2016/02/09 : non à la ligne suivante : il faut que METASS soit
*     égale à 5 (voir remarque dans makpr2)
*      METASS=4
      METASS=5
*
      CALL KRES3(MATRIK,MATASS,MRENU,MMULAG,METASS,
     $     0,0,0,0,0,
     $     0,.FALSE.,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C
C Extraction de la diagonale
C
      SEGACT MATRIK
      AMORS=MATRIK.KIDMAT(4)
      AISA=MATRIK.KIDMAT(5)
*
      SEGACT AMORS
      SEGACT AISA
      NTTDDL=AMORS.IA(/1)-1
      NBVA=NTTDDL
      SEGINI IDIAGO
      IF (IOPT.EQ.2) THEN
         SEGINI ISCAR
      ENDIF
      DO ITTDDL=1,NTTDDL
         JSTRT=AMORS.IA(ITTDDL)
         JSTOP=AMORS.IA(ITTDDL+1)-1
         DO J=JSTRT,JSTOP
            JTTDDL=AMORS.JA(J)
            IF (JTTDDL.EQ.ITTDDL) THEN
               IDIAGO.A(ITTDDL)=AISA.A(J)
            ENDIF
            IF (IOPT.EQ.2) THEN
               VAL=AISA.A(J)
               ISCAR.A(ITTDDL)=ISCAR.A(ITTDDL)+(VAL*VAL)
            ENDIF
         ENDDO
      ENDDO
      SEGSUP AISA
      SEGSUP AMORS
*
      IF (IOPT.EQ.2) THEN
         DO ITTDDL=1,NTTDDL
            VAL=ISCAR.A(ITTDDL)
            IF (VAL.LE.SQRT(XPETIT)) THEN
               WRITE(IOIMP,*) 'La ligne ',ITTDDL,
     $              ' de la matrice est nulle : ', VAL
               GOTO 9999
            ENDIF
            IDIAGO.A(ITTDDL)=IDIAGO.A(ITTDDL)/VAL
         ENDDO
      ENDIF
      IF (IOPT.EQ.2) THEN
         SEGSUP ISCAR
      ENDIF
C
C Transformation en chpoint
C
      CALL XDISP(MATRIK,IDIAGO,MCHDIA,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      SEGSUP IDIAGO
      SEGSUP MATRIK
      CALL ECROBJ('CHPOINT ',MCHDIA)
*
* Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in exdiag.eso'
*  153 2
* Opération illicite dans ce contexte
      CALL ERREUR(153)
      RETURN
*
* End of EXDIAG
*
      END