Sources Esope | Cast3M

kres4
C KRES4     SOURCE    GOUNAND   22/08/25    21:15:06     11434          
      SUBROUTINE KRES4(MATRIK,KCLIM,KSMBR,
     $     ISCAL,
     $     MCHSOL,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : KRES4
C DESCRIPTION : Résolution d'un système par une méthode directe
C               (Factorisation LDU sans pivoting).
C
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : VERMAT, MESMBR, MELIM, CLMORS, TRIALU, REZOLU,
C                    XDISP
C APPELES (E/S)    : INFMAT
C APPELE PAR       : KRES2
C***********************************************************************
C ENTREES            : MATRIK, KCLIM, KSMBR
C ENTREES/SORTIES    : -
C SORTIES            : MCHSOL
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 14/04/2000, version initiale
C HISTORIQUE : v1, 14/04/2000, création
C HISTORIQUE : 06/04/04 : Scaling
C HISTORIQUE :
C HISTORIQUE :
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 SMCHPOI
      POINTEUR KCLIM.MCHPOI
      POINTEUR KSMBR.MCHPOI
      POINTEUR MCHSOL.MCHPOI
-INC SMMATRIK
      POINTEUR INCX.IZA
      POINTEUR KS2B.IZA
      POINTEUR KMORS.PMORS
      POINTEUR KISA.IZA
      POINTEUR AMORS.PMORS
      POINTEUR AISA.IZA
      POINTEUR NORMP.IZA
      POINTEUR NORMD.IZA
*
      INTEGER IMPR,IRET
*
* Executable statements
*
      IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans kres4.eso'
C
C On vérifie que la matrice est correctement assemblée
C
      CALL VERMAT(MATRIK,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C Pas besoin d'estimation de départ pour obtenir la solution
      INCX=0
C
C On transforme le chpoint second membre en vecteur second membre
C In MESMBR : SEGINI KS2B
      CALL MESMBR(MATRIK,KSMBR,
     $     KS2B,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C
C On applique les conditions aux limites
C
C In MELIM : SEGINI AMORS
C            SEGINI AISA
      CALL MELIM(MATRIK,KCLIM,
     $     INCX,KS2B,
     $     AMORS,AISA,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C Eventuellement, on élimine les 0.d0 de la matrice
C On ne le fait que pour les méthodes directes (où le gain en temps est
C important) .
C Il faudrait régler les problèmes de précision pour les termes
C proches de zéro pouvant engendrer des pertes de symétrie
C de la matrice...
C Il y a aussi des problèmes lorsqu'un
C terme vaut zéro à une itération et autre chose aux suivantes.
C L'idéal serait de ne virer que les 0.D0 issus des conditions
C aux limites...
      CALL CLMORS(AMORS,AISA,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      IF (ISCAL.GT.0) THEN
C
C Calcul des normes primales (colonnes) et duales (lignes)
C de la matrice. Norme = norme L2, soit :
C {\sum_{i ou j} a_{ij}^2}^{1/2}
C
         CALL NORMAT(AMORS,AISA,ISCAL,NORMP,NORMD,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
C
C On norme la matrice : attention modification...
C
         CALL NORMAM(AMORS,AISA,NORMP,NORMD,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
C
C On norme le second membre : attention modification...
C
         CALL NORMV2(KS2B,NORMD,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
      ENDIF
C
C On donne des infos sur la matrice
C
      CALL INFMAT(AMORS,AISA,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C
C Factorisation LDU de la matrice
C
      CALL TRIALU(MATRIK,AMORS,AISA,
     $     IDMAT,
     $     IMPR,IRET)
      if (ierr.ne.0) return
      IF (IRET.NE.0) GOTO 9999
C
C Obtention de la solution (montée-descente)
C
      CALL REZOLU(IDMAT,
     $     KS2B,
     $     IMPR,IRET)
      if (ierr.ne.0) return
      IF (IRET.NE.0) GOTO 9999
      INCX=KS2B
      IF (ISCAL.GT.0) THEN
C
C On dénorme le vecteur solution : attention modification...
C
         CALL NORMV2(INCX,NORMP,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
      ENDIF
C
C Transformation du vecteur-solution en chpoint
C
      CALL XDISP(MATRIK,INCX,MCHSOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C
C Suppression des objets temporaires
C
      SEGSUP,AMORS
      SEGSUP,AISA
      SEGSUP INCX
*
* Normal termination
*
      IRET=0
      RETURN
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in kres4.eso'
      RETURN
*
* End of KRES4
*
      END