Esope source files

   1 : C KRES4     SOURCE    MAGN      09/07/30    21:18:52     6458           
   2 :       SUBROUTINE KRES4(MATRIK,KCLIM,KSMBR,
   3 :      $     ISCAL,
   4 :      $     MCHSOL,
   5 :      $     IMPR,IRET)
   6 :       IMPLICIT REAL*8 (A-H,O-Z)
   7 :       IMPLICIT INTEGER (I-N)
   8 : C***********************************************************************
   9 : C NOM         : KRES4
  10 : C DESCRIPTION : Résolution d'un système par une méthode directe
  11 : C               (Factorisation LDU sans pivoting).
  12 : C               
  13 : C
  14 : C LANGAGE     : ESOPE
  15 : C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  16 : C               mél : gounand@semt2.smts.cea.fr
  17 : C***********************************************************************
  18 : C APPELES          : VERMAT, MESMBR, MELIM, CLMORS, TRIALU, REZOLU,
  19 : C                    XDISP
  20 : C APPELES (E/S)    : INFMAT
  21 : C APPELE PAR       : KRES2
  22 : C***********************************************************************
  23 : C ENTREES            : MATRIK, KCLIM, KSMBR
  24 : C ENTREES/SORTIES    : -
  25 : C SORTIES            : MCHSOL
  26 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  27 : C***********************************************************************
  28 : C VERSION    : v1, 14/04/2000, version initiale
  29 : C HISTORIQUE : v1, 14/04/2000, création
  30 : C HISTORIQUE : 06/04/04 : Scaling
  31 : C HISTORIQUE :
  32 : C HISTORIQUE :
  33 : C***********************************************************************
  34 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
  35 : C en cas de modification de ce sous-programme afin de faciliter
  36 : C la maintenance !
  37 : C***********************************************************************
  38 : -INC CCOPTIO
  39 : -INC SMCHPOI
  40 :       POINTEUR KCLIM.MCHPOI
  41 :       POINTEUR KSMBR.MCHPOI
  42 :       POINTEUR MCHSOL.MCHPOI
  43 : -INC SMMATRIK
  44 :       POINTEUR INCX.IZA
  45 :       POINTEUR KS2B.IZA
  46 :       POINTEUR KMORS.PMORS
  47 :       POINTEUR KISA.IZA
  48 :       POINTEUR AMORS.PMORS
  49 :       POINTEUR AISA.IZA
  50 :       POINTEUR NORMP.IZA
  51 :       POINTEUR NORMD.IZA
  52 : *
  53 :       INTEGER IMPR,IRET
  54 : *     
  55 : * Executable statements
  56 : *
  57 :       IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans kres4.eso' 
  58 : C
  59 : C On vérifie que la matrice est correctement assemblée
  60 : C
  60 :       CALL VERMAT(MATRIK,IMPR,IRET)
  61 :       IF (IRET.NE.0) GOTO 9999
  62 : C Pas besoin d'estimation de départ pour obtenir la solution
  63 :       INCX=0
  64 : C
  65 : C On transforme le chpoint second membre en vecteur second membre
  66 : C In MESMBR : SEGINI KS2B
  66 :       CALL MESMBR(MATRIK,KSMBR,
  67 :      $     KS2B,
  68 :      $     IMPR,IRET)
  69 :       IF (IRET.NE.0) GOTO 9999
  70 : C
  71 : C On applique les conditions aux limites
  72 : C
  73 : C In MELIM : SEGINI AMORS 
  74 : C            SEGINI AISA
  74 :       CALL MELIM(MATRIK,KCLIM,
  75 :      $     INCX,KS2B,
  76 :      $     AMORS,AISA,
  77 :      $     IMPR,IRET)
  78 :       IF (IRET.NE.0) GOTO 9999
  79 : C Eventuellement, on élimine les 0.d0 de la matrice
  80 : C On ne le fait que pour les méthodes directes (où le gain en temps est
  81 : C important) .
  82 : C Il faudrait régler les problèmes de précision pour les termes
  83 : C proches de zéro pouvant engendrer des pertes de symétrie
  84 : C de la matrice...
  85 : C Il y a aussi des problèmes lorsqu'un
  86 : C terme vaut zéro à une itération et autre chose aux suivantes.
  87 : C L'idéal serait de ne virer que les 0.D0 issus des conditions
  88 : C aux limites...
  88 :       CALL CLMORS(AMORS,AISA,IMPR,IRET)
  89 :       IF (IRET.NE.0) GOTO 9999
  90 :       IF (ISCAL.GT.0) THEN
  91 : C
  92 : C Calcul des normes primales (colonnes) et duales (lignes) 
  93 : C de la matrice. Norme = norme L2, soit :
  94 : C {\sum_{i ou j} a_{ij}^2}^{1/2}
  95 : C
  95 :          CALL NORMAT(AMORS,AISA,ISCAL,NORMP,NORMD,IMPR,IRET)
  96 :          IF (IRET.NE.0) GOTO 9999
  97 : C
  98 : C On norme la matrice : attention modification...
  99 : C
  99 :          CALL NORMAM(AMORS,AISA,NORMP,NORMD,IMPR,IRET)
 100 :          IF (IRET.NE.0) GOTO 9999
 101 : C
 102 : C On norme le second membre : attention modification...
 103 : C
 103 :          CALL NORMV2(KS2B,NORMD,IMPR,IRET)
 104 :          IF (IRET.NE.0) GOTO 9999
 105 :       ENDIF
 106 : C
 107 : C On donne des infos sur la matrice
 108 : C
 108 :       CALL INFMAT(AMORS,AISA,IMPR,IRET)
 109 :       IF (IRET.NE.0) GOTO 9999
 110 : C
 111 : C Factorisation LDU de la matrice
 112 : C
 112 :       CALL TRIALU(MATRIK,AMORS,AISA,
 113 :      $     IDMAT,
 114 :      $     IMPR,IRET)
 115 :       IF (IRET.NE.0) GOTO 9999
 116 : C
 117 : C Obtention de la solution (montée-descente)
 118 : C
 118 :       CALL REZOLU(IDMAT,
 119 :      $     KS2B,
 120 :      $     IMPR,IRET)
 121 :       IF (IRET.NE.0) GOTO 9999
 122 :       INCX=KS2B
 123 :       IF (ISCAL.GT.0) THEN
 124 : C
 125 : C On dénorme le vecteur solution : attention modification...
 126 : C
 126 :          CALL NORMV2(INCX,NORMP,IMPR,IRET)
 127 :          IF (IRET.NE.0) GOTO 9999
 128 :       ENDIF
 129 : C
 130 : C Transformation du vecteur-solution en chpoint
 131 : C
 131 :       CALL XDISP(MATRIK,INCX,MCHSOL,IMPR,IRET)
 132 :       IF (IRET.NE.0) GOTO 9999
 133 : C
 134 : C Suppression des objets temporaires
 135 : C
 136 :       SEGSUP,AMORS
 137 :       SEGSUP,AISA
 138 :       SEGSUP INCX
 139 : *     
 140 : * Normal termination
 141 : *     
 142 :       IRET=0
 143 :       RETURN
 144 : *     
 145 : * Error handling
 146 : *     
 147 :  9999 CONTINUE
 148 :       IRET=1
 149 :       WRITE(IOIMP,*) 'An error was detected in kres4.eso'
 150 :       RETURN
 151 : *     
 152 : * End of KRES4
 153 : *     
 154 :       END
 155 : 
 156 : 
 157 : 
 158 : 
 159 : 
 160 : 
 161 : 
 162 :  
 163 :  
 164 :  
 165 :  
 166 :  
 167 :  
 168 :  
 169 :  
 170 :