Télécharger la plaquette Cast3M


fr en

Télécharger kres8.eso

Retour à la liste

Numérotation des lignes : 

kres8

  1. C KRES8     SOURCE    GOUNAND   25/04/30    21:15:17     12258          
  2.       SUBROUTINE KRES8(MRIGID,KSMBR,INORMU,
  3.      $     KTYPI,ITER,RESID,ICALRS,IRSTRT,LBCG,BRTOL,IDDOT,IMVEC,
  4.      $     IORINC,MLAG1,MLAG2,
  5.      $     KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
  6.      $     KTIME,LTIME,LDUMP,ISMOOT,
  7.      $     MCHSOL,LRES,LNMV,ICVG,IMPR)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9.       IMPLICIT INTEGER (I-N)
  10. C***********************************************************************
  11. C NOM         : KRES8
  12. C DESCRIPTION : - Assemblage par RESOU
  13. C               - Conversion au format Morse de la matrice
  14. C               - Conversion du second membre en MVECTD
  15. C               - Numerotation Bloc ou multiplicateur de Lagrange
  16. C                  (KRES13-KRES14)
  17. C               - Construction du préconditionneur
  18. C               - Appel des solveurs itératifs
  19. C               - Conversion du résultat en CHPOINT
  20. C
  21. C
  22. C LANGAGE     : ESOPE
  23. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  24. C               mél : gounand@semt2.smts.cea.fr
  25. C***********************************************************************
  26. C VERSION    : v1, 04/08/2011, version initiale
  27. C HISTORIQUE : v1, 04/08/2011, création
  28. C HISTORIQUE :
  29. C HISTORIQUE :
  30. C***********************************************************************
  31.  
  32. -INC PPARAM
  33. -INC CCOPTIO
  34. -INC CCREEL
  35. C On inclue SMCOORD car MCHSOL doit avoir la configuration courante
  36. -INC SMCOORD
  37. -INC SMCHPOI
  38.       POINTEUR MCHSOL.MCHPOI
  39. -INC SMRIGID
  40. -INC SMVECTD
  41.       POINTEUR ASMBR.MVECTD
  42.       POINTEUR ISMBR.MVECTD
  43.       POINTEUR INCX.MVECTD
  44.       POINTEUR IR.MVECTD
  45. -INC SMMATRI
  46.       SEGMENT PMORS
  47.       INTEGER IA (NTT+1)
  48.       INTEGER JA (NJA)
  49.       ENDSEGMENT
  50.       POINTEUR PMS1.PMORS,PMS2.PMORS
  51.       POINTEUR KMORS.PMORS,AMORS.PMORS
  52. C Segment de stokage
  53.       SEGMENT IZA
  54.       REAL*8  A(NBVA)
  55.       ENDSEGMENT
  56.       POINTEUR IZA1.IZA,IZA2.IZA,IZAU.IZA,IZAL.IZA,ISA.IZA
  57.       POINTEUR KIZA.IZA,AIZA.IZA
  58.  
  59. -INC SMLENTI
  60.       POINTEUR KTYP.MLENTI
  61.       POINTEUR KNOD.MLENTI
  62.       POINTEUR KRINC.MLENTI
  63.       POINTEUR IWORK.MLENTI
  64.       POINTEUR LAGRAN.MLENTI
  65.       POINTEUR JORDRE.MLENTI,JORTMP.MLENTI
  66.       POINTEUR IPERM.MLENTI,JPETMP.MLENTI
  67.       POINTEUR JPERM.MLENTI
  68.       POINTEUR IPBLOC.MLENTI
  69. -INC SMLMOTS
  70.       POINTEUR IORINC.MLMOTS
  71.       POINTEUR IORINU.MLMOTS
  72. -INC SMTABLE
  73.       POINTEUR KTIME.MTABLE
  74.       DIMENSION ITTIME(4)
  75.       CHARACTER*(LOCHPO) CHCOMP
  76.       CHARACTER*16 CHARI
  77.       CHARACTER*1 CCOMP
  78.       LOGICAL LTIME,LOGII,LMG,LDUMP,LVERIF
  79.       REAL*8 GNRM2
  80. C     ..
  81. C     .. External subroutines and functions..
  82. *inutile      EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
  83.  
  84.       LVERIF=.FALSE.
  85.       IVALI=0
  86.       XVALI=0.D0
  87.       LOGII=.FALSE.
  88.       IRETI=0
  89.       XVALR=0.D0
  90.       IRETR=0
  91. C
  92. C Executable statements
  93. C
  94.       IF (LTIME) THEN
  95.          CALL CRTABL(KTIME)
  96.          call timespv(ittime,oothrd)
  97.          ITI1=(ITTIME(1)+ITTIME(2))/10
  98.       ELSE
  99.          KTIME=0
  100.       ENDIF
  101. C
  102. C  CAS PARTICULIER : Si la matrice est vide (toutes les inconnues
  103. C      éliminées, par exemple)
  104. C
  105.       SEGACT MRIGID
  106.       IF (IRIGEL(/2).EQ.0) THEN
  107.          NSOUPO=0
  108.          NAT=0
  109.          SEGINI MCHSOL
  110.          SEGDES MCHSOL
  111.          ICVG=0
  112.          LNMV=0
  113.          LRES=0
  114.          IF (LTIME) THEN
  115.             call timespv(ittime,oothrd)
  116.             ITI2=(ITTIME(1)+ITTIME(2))/10
  117.             CHARI='MATVIDE'
  118.             CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  119.      $           'ENTIER  ',ITI2-ITI1,XVALR,CHARR,LOGIR,IRETR)
  120.             SEGDES KTIME
  121.          ENDIF
  122.          SEGDES MRIGID
  123.          RETURN
  124.       ENDIF
  125. C
  126. C               - Assemblage par RESOU
  127. C
  128. C  old INORMU=1 : Normalisation des mutiplicateurs de Lagrange
  129. *     INORMU est transmis à la subroutine
  130. *     Le problème est que si MRIGID est deja assemblée, INORMU n'est pas
  131. *     pris en compte... mais où le stocker ??
  132.       CALL KRES9(MRIGID,INORMU)
  133.       IF (IERR.NE.0) RETURN
  134.       IF (LTIME) THEN
  135.          call timespv(ittime,oothrd)
  136.          ITI2=(ITTIME(1)+ITTIME(2))/10
  137.       ENDIF
  138. C
  139. C               - Conversion au format Morse de la matrice
  140. C
  141.       CALL KRES10(MRIGID,KMORS,KIZA)
  142.       IF (IERR.NE.0) RETURN
  143.       IF (LTIME) THEN
  144.          call timespv(ittime,oothrd)
  145.          ITI3=(ITTIME(1)+ITTIME(2))/10
  146.       ENDIF
  147. C
  148. C On donne des infos sur la matrice
  149. C
  150. *      SEGACT MRIGID
  151. *     ICHOLX=ICHOLE
  152. **    INFDDL.ESO est dans ~/triou/p1nc
  153. **      CALL INFDDL(ICHOLX)
  154. C      WRITE(IOIMP,*) 'IMPR=',IMPR
  155.       CALL INFMAT(KMORS,KIZA,IMPR,IRET)
  156. C      IF (IRET.NE.0) GOTO 9999
  157. C      WRITE(IOIMP,*) 'Apres KRES10'
  158. C      WRITE(IOIMP,*) 'KMORS=',KMORS
  159. C      WRITE(IOIMP,*) 'KIZA=',KIZA
  160.  
  161. C
  162. C               - Conversion du second membre en MVECTD
  163. C                 et initialisation du résultat
  164. C
  165.       SEGACT MRIGID
  166.       ICHOLX=ICHOLE
  167.       ISECO=KSMBR
  168. C On ne vérifie pas que le second membre doit être dans le dual
  169.       NOID=1
  170.       CALL CHVNS(ICHOLX,ISECO,ISMBR,NOID)
  171.       IF (IERR.NE.0) RETURN
  172.       IF (LTIME) THEN
  173.          call timespv(ittime,oothrd)
  174.          ITI4=(ITTIME(1)+ITTIME(2))/10
  175.       ENDIF
  176.  
  177. C      SEGACT ISMBR
  178. C      WRITE(IOIMP,*) 'Second membre sous forme vecteur'
  179. C      INC=ISMBR.VECTBB(/1)
  180. C      WRITE(IOIMP,*) '     ISMBR, INC=',INC
  181. C      WRITE(IOIMP,2022) (ISMBR.VECTBB(II),II=1,ISMBR.VECTBB(/1))
  182. C
  183. C     Gestion normalisation Lagrange (repris de MONDES)
  184. C
  185. *      IF (INORMU.EQ.1) THEN
  186.       SEGACT ISMBR*MOD
  187.       MMATRI=ICHOLE
  188.       SEGACT MMATRI
  189.       IF(IDNORD.GT.0) THEN
  190.          MDNO1=IDNORD
  191.       ELSE
  192.          MDNO1=IDNORM
  193.       ENDIF
  194.       SEGACT MDNO1
  195.       INC=MDNO1.DNOR(/1)
  196.       DO 45 I=1,INC
  197.          ISMBR.VECTBB(I)=ISMBR.VECTBB(I)*MDNO1.DNOR(I)
  198.  45   CONTINUE
  199. *      ENDIF
  200. C
  201. C     Placer correctement les multiplicateurs de Lagrange
  202. C     ou réordonnancer par blocs pour le multigrille
  203. C
  204.       LMG=(KTYPI.EQ.7.OR.KTYPI.EQ.8.OR.KTYPI.EQ.10.OR.KTYPI.EQ.11)
  205. *      write(ioimp,*) 'LMG=',LMG
  206.       IF (LMG) THEN
  207.          CALL KRES13(MRIGID,KMORS,KIZA,ISMBR,IORINC,KTYPI,
  208.      $        IPERM,JPERM,IPBLOC)
  209.          IF (IERR.NE.0) RETURN
  210.       ELSE
  211.          CALL KRES14(MRIGID,KMORS,KIZA,ISMBR,MLAG1,MLAG2,
  212.      $        IPERM,JPERM)
  213.          IF (IERR.NE.0) RETURN
  214.          IPBLOC=0
  215.       ENDIF
  216. *     segprt,iperm
  217. *     segprt,jperm
  218.       if (iperm.ne.0) then
  219.          SEGINI,AMORS=KMORS
  220.          SEGINI,AIZA=KIZA
  221.          SEGINI,ASMBR=ISMBR
  222. *      write(ioimp,*) 'Avant permutation'
  223. *      CALL ECMORS(AMORS,AIZA,3)
  224.          job=1
  225.          nbddl=AMORS.IA(/1)-1
  226.          CALL DPERM(nbddl,kiza.a,kmors.ja,kmors.ia,
  227.      $        aiza.a,amors.ja,amors.ia,iperm.lect,jperm.lect,job)
  228.          CALL DVPERM(nbddl,asmbr.vectbb,iperm.lect)
  229. *      write(ioimp,*) 'Apres permutation'
  230. *      CALL ECMORS(AMORS,AIZA,3)
  231.          segsup iperm
  232.       else
  233.          AMORS=KMORS
  234.          AIZA=KIZA
  235.          ASMBR=ISMBR
  236.       endif
  237. *
  238. *     ordonnancement des colonnes
  239. *     En multigrille, on n'est pas oblige mais c'est pas tres long de
  240. *     toutes facons
  241. *      if (.NOT.lmg) then
  242.          CALL KRES15(AMORS,AIZA)
  243. *      endif
  244. *      write(ioimp,*) 'Apres ordonnancement colonnes'
  245. *      CALL ECMORS(AMORS,AIZA,3)
  246.       IF (IERR.NE.0) RETURN
  247. C
  248. C     - Construction du préconditionneur (repris sur kres5)
  249. C     - Appel des solveurs itératifs
  250. C
  251. C         CALL ECMORS(AMORS,AIZA,4)
  252. C      SEGACT ASMBR
  253. C      WRITE(IOIMP,*) 'Second membre sous forme vecteur'
  254. C      INC=ASMBR.VECTBB(/1)
  255. C      WRITE(IOIMP,*) '     ASMBR, INC=',INC
  256. C      WRITE(IOIMP,2022) (ASMBR.VECTBB(II),II=1,ASMBR.VECTBB(/1))
  257. C Solveur Direct
  258.       IF (KTYPI.EQ.1) THEN
  259.          SEGINI,INCX=ASMBR
  260.          CALL KRES12(AMORS,AIZA,INCX,
  261. C         CALL KRES12(KMORS,KIZA,ISMBR,
  262.      $        KTIME,LTIME,
  263.      $        INCX,LRES,LNMV,ICVG,IMPR,INODET)
  264.       ELSE
  265. C Solveur Itératif
  266.          CALL KRES11(AMORS,AIZA,ASMBR,IPBLOC,
  267.      $        KTYPI,ITER,RESID,ICALRS,IRSTRT,LBCG,BRTOL,IDDOT,IMVEC,
  268.      $        KPREC,RXMILU,RXILUP,XLFIL,XDTOL,XSPIV,
  269.      $        KTIME,LTIME,LDUMP,ISMOOT,
  270.      $        INCX,LRES,LNMV,ICVG,IMPR,INODET)
  271. C      WRITE(IOIMP,*) 'Apres KRES11'
  272.       ENDIF
  273.       IF(IERR.NE.0) RETURN
  274. *
  275.       IF (LVERIF) THEN
  276. * Verif résidu
  277.          SEGINI,IR=ASMBR
  278.          SEGACT INCX
  279.          SEGACT AMORS
  280.          SEGACT AIZA
  281. C                    r(0)=b-Ax
  282.          CALL GMOMV(IMVEC,'N',-1.D0,AMORS,AIZA,INCX,1.D0,IR)
  283.          RNRM2 = GNRM2(IR)
  284.          WRITE(IOIMP,*) '||R||_1 =',RNRM2
  285.          SEGSUP,IR
  286.       ENDIF
  287. *
  288.  
  289.  
  290. *
  291.       SEGACT INCX*MOD
  292.       IF (JPERM.NE.0) THEN
  293.          CALL DVPERM(nbddl,incx.vectbb,jperm.lect)
  294.          segsup jperm
  295.          segsup amors
  296.          segsup aiza
  297.          segsup asmbr
  298.       ENDIF
  299.       if (ipbloc.ne.0) segsup ipbloc
  300. C      SEGACT INCX
  301. C      WRITE(IOIMP,*) 'Inconnue sous forme vecteur'
  302. C      INC=INCX.VECTBB(/1)
  303. C      WRITE(IOIMP,*) '     INCX, INC=',INC
  304. C      WRITE(IOIMP,2022) (INCX.VECTBB(II),II=1,INCX.VECTBB(/1))
  305. C      IF(IERR.NE.0) RETURN
  306.       IF (LVERIF) THEN
  307. C                    r(0)=b
  308.          SEGINI,IR=ISMBR
  309.          SEGACT KMORS
  310.          SEGACT KIZA
  311. C                    r(0)=b-Ax
  312.          CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,IR)
  313.          RNRM2 = GNRM2(IR)
  314.          WRITE(IOIMP,*) '||R||_2=',RNRM2
  315.       ENDIF
  316. C
  317. C   Pour l'instant, on ne stocke pas la matrice Morse
  318. C
  319.       SEGSUP KMORS
  320.       SEGSUP KIZA
  321. C
  322. C     Gestion normalisation Lagrange (repris de MONDES)
  323. C     + égalité multiplicateurs
  324. *      IF (INORMU.EQ.1) THEN
  325.          MMATRI=ICHOLE
  326.          SEGACT MMATRI
  327.          MDNOR=IDNORM
  328.          SEGACT MDNOR
  329.          INC=DNOR(/1)
  330.          DO 35 I=1,INC
  331.             INCX.VECTBB(I)=INCX.VECTBB(I)*DNOR(I)
  332.  35      CONTINUE
  333.          SEGDES MDNOR
  334.          MILIGN=IILIGN
  335.          SEGACT,MILIGN
  336.          DO 36 I = 1, INC
  337.             if (ITTR(I).ne.0) then
  338. *            write (6,*) ' dans mondes ',i,ittr(i)
  339.                if (incx.vectbb(i).eq.0.d0
  340.      $              .or.incx.vectbb(ittr(i)).eq.0.d0) then
  341. *      write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k)
  342.                   incx.vectbb(i)=0.d0
  343.                   incx.vectbb(ittr(i))=0.d0
  344.                   goto 36
  345.                endif
  346.                incx.vectbb(i)=(incx.vectbb(i)+incx.vectbb(ittr(i)))/2
  347.                incx.vectbb(ittr(i))=incx.vectbb(i)
  348.             endif
  349.  36      CONTINUE
  350.          SEGDES MILIGN
  351.          SEGDES MMATRI
  352.          SEGDES INCX
  353. *      ENDIF
  354. C
  355. C
  356. C      SEGACT INCX
  357. C      WRITE(IOIMP,*) 'Inconnue sous forme vecteur'
  358. C      INC=INCX.VECTBB(/1)
  359. C      WRITE(IOIMP,*) '     INCX, INC=',INC
  360. C      WRITE(IOIMP,2022) (INCX.VECTBB(II),II=1,INCX.VECTBB(/1))
  361. C      IF(IERR.NE.0) RETURN
  362.  
  363.       IF (LTIME) THEN
  364.          call timespv(ittime,oothrd)
  365.          ITI5=(ITTIME(1)+ITTIME(2))/10
  366.       ENDIF
  367. C
  368. C               - Conversion du résultat en CHPOINT
  369. C
  370.       CALL VCH1(ICHOLX,INCX,MCHSOL,MRIGID)
  371. C      WRITE(IOIMP,*) 'Apres VCH1'
  372.       IF(IERR.NE.0) RETURN
  373.       IF (LTIME) THEN
  374.          call timespv(ittime,oothrd)
  375.          ITI6=(ITTIME(1)+ITTIME(2))/10
  376.          CHARI='RESO ASS+RENU'
  377.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  378.      $                     'ENTIER  ',ITI2-ITI1,XVALR,CHARR,LOGIR,IRETR)
  379.          CHARI='CONVMORS'
  380.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  381.      $                     'ENTIER  ',ITI3-ITI2,XVALR,CHARR,LOGIR,IRETR)
  382. C         CHARI='CONVSMB '
  383. C         CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  384. C     $                     'ENTIER  ',ITI4-ITI3,XVALR,CHARR,LOGIR,IRETR)
  385.          IF (KTYPI.EQ.1) THEN
  386.             CHARI='KRES FAC+RESO'
  387.          ELSE
  388.             CHARI='KRES PRE+RESO'
  389.          ENDIF
  390.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  391.      $                     'ENTIER  ',ITI5-ITI4,XVALR,CHARR,LOGIR,IRETR)
  392. C         CHARI='CONVINC'
  393. C         CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  394. C     $                     'ENTIER  ',ITI6-ITI5,XVALR,CHARR,LOGIR,IRETR)
  395.          CHARI='TOTAL   '
  396.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  397.      $                     'ENTIER  ',ITI6-ITI1,XVALR,CHARR,LOGIR,IRETR)
  398.          SEGDES KTIME
  399.       ENDIF
  400. C     Le solveur direct surcharge le second membre
  401.       IF (ISMBR.NE.INCX) SEGSUP ISMBR
  402.       SEGSUP INCX
  403.       SEGDES MRIGID
  404. C
  405. C Normal termination
  406. C
  407.       RETURN
  408. C
  409. C Error Handling
  410. C
  411.  9999 CONTINUE
  412.       MOTERR(1:8)='KRES8   '
  413.       CALL ERREUR(1127)
  414.       RETURN
  415. C
  416. C Format handling
  417. C
  418.  2022 FORMAT(10(1X,1PG12.5))
  419. C
  420. C End of subroutine KRES8
  421. C
  422.       END
  423.  
  424.  

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