Télécharger la plaquette Cast3M


fr en

Télécharger krig.eso

Retour à la liste

Numérotation des lignes : 

krig

  1. C KRIG      SOURCE    FD218221  25/09/03    21:15:02     12354          
  2.       SUBROUTINE KRIG
  3.  
  4. C     Typages implicites habituels
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.  
  8. C     Les includes necessaires
  9. -INC PPARAM
  10. -INC CCOPTIO
  11. -INC CCREEL
  12. -INC SMLMOTS
  13. -INC SMLREEL
  14. -INC SMEVOLL
  15. -INC SMELEME
  16. -INC SMCHPOI
  17. -INC SMTABLE
  18.  
  19. C     Un segment pour stocker les donnees
  20. C     NC : nombre de coordonnees des points
  21. C     MM : nombre de points de mesure
  22. C     NN : nombre de points cibles
  23. C     MCOO(i,k) : k-eme coordonnee du i-eme point de mesure
  24. C     MVAL(i,1) : valeur fonction au  i-eme point de mesure
  25. C     CCOO(j,k) : k-eme coordonnee du j-eme point cible
  26. C     DM(i1,i2) : distance entre les points de mesure i1 et i2
  27. C     DC(i,j)   : distance entre le point de mesure i et le point cible j
  28. C     VM(i1,i2) : variance (ou covariance) entre les points de mesure i1 et i2
  29. C                 la derniere ligne est (1 1 1 ... 1 0) pour assurer le krigeage ordinaire
  30. C     VC(i,j)   : Avant appel a DGESV : variance (ou covariance) entre le points mesure i
  31. C                                       et point cible j la derniere ligne est (1 1 1 ... 1)
  32. C                                       pour assurer le krigeage ordinaire
  33. C                 Apres appel a DGESV : poids d'interpolation pour la i-eme mesure
  34. C                                       pour interpoler au j-eme point cible
  35. C                                       la derniere contient les multiplicateurs de Lagrange
  36. C                                       pour assurer le krigeage ordinaire
  37.       SEGMENT SKRIG
  38.           REAL*8 MCOO(MM,NC)
  39.           REAL*8 MVAL(MM,1)
  40.           REAL*8 CCOO(NN,NC)
  41.           REAL*8 DM(MM,MM)
  42.           REAL*8 DC(MM,NN)
  43.           REAL*8 VM(MM+1,MM+1)
  44.           REAL*8 VC(MM+1,NN)
  45.           REAL*8 WC(MM+1,NN)
  46.           INTEGER IPIV(MM+1)
  47.       ENDSEGMENT
  48.  
  49. C     Quelques SEGMENTs LISTREEL
  50.       POINTEUR MLREE4.MLREEL
  51.  
  52. C     Quelques objets utiles
  53.       PARAMETER(XUN=1.D0)
  54.       CHARACTER*4 MOT1,COMP
  55.       CHARACTER*8 TYPOBJ
  56.       LOGICAL BVALE
  57.  
  58.  
  59.  
  60. C----------------------------------------------------------------------C
  61. C             ACQUISITION DES DONNEES D'ENTREE UTILISATEUR             C
  62. C----------------------------------------------------------------------C
  63.  
  64. C     La table (ITAB)
  65.       CALL LIROBJ('TABLE',ITAB,1,IRETOU)
  66.       IF(IERR.NE.0) RETURN
  67.  
  68. C     Noms des dimensions (LISTMOTS : MLMOT1)
  69.       CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'COORDONNEES',BIND,IOIND,
  70.      &            'LISTMOTS',IVALE,XVALE,' ',BVALE,MLMOT1)
  71.       IF (IERR.NE.0) GOTO 999
  72.       SEGACT MLMOT1
  73.  
  74. C     Mesures (CHPOINT : MCHPO1)
  75.       CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'MESURES',BIND,IOIND,
  76.      &            'CHPOINT',IVALE,XVALE,' ',BVALE,MCHPO1)
  77.       IF (IERR.NE.0) GOTO 999
  78.       SEGACT MCHPO1
  79. C     Maillage de POI1 support des mesures (IPT1)
  80.       CALL EXTR21(MCHPO1,1,IPT1)
  81. C     Nombre de points de mesure (M)
  82.       M=IPT1.NUM(/2)
  83.  
  84. C     Nom de la composante a kirger (MOT : COMP)
  85.       CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'COMPOSANTE',BIND,IOIND,
  86.      &            'MOT     ',IVALE,XVALE,COMP,BVALE,IOVALE)
  87.       IF (IERR.NE.0) GOTO 999
  88.  
  89. C     Cibles (CHPOINT : MCHPO2)
  90.       CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'CIBLES',BIND,IOIND,
  91.      &            'CHPOINT ',IVALE,XVALE,' ',BVALE,MCHPO2)
  92.       IF (IERR.NE.0) GOTO 999
  93.       SEGACT MCHPO2
  94. C     Maillage de POI1 support des cibles (IPT2)
  95.       CALL EXTR21(MCHPO2,1,IPT2)
  96. C     Nombre de points cibles (N)
  97.       N=IPT2.NUM(/2)
  98.  
  99. C     Variogramme ? (IVAR=1) ou covariogramme ? (IVAR=2) (EVOLUTIO : MEVOL1)
  100.       IVAR=0
  101.       TYPOBJ='        '
  102.       CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'VARIOGRAMME',BIND,IOIND,
  103.      &            TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  104.       IF (IERR.NE.0) GOTO 999
  105.       IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  106.           IVAR=1
  107.       ELSE
  108.           CALL ACCTAB(ITAB,'MOT     ',IIND,XIND,'COVARIOGRAMME',BIND,
  109.      &          IOIND,TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  110.           IF (IERR.NE.0) GOTO 999
  111.           IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  112.               IVAR=2
  113.           ELSE
  114.               MOTERR(1:40)='VARIOGRAMME ou COVARIOGRAMME            '
  115.               CALL ERREUR(535)
  116.               RETURN
  117.           ENDIF
  118.       ENDIF
  119. C     Recuperation des 2 LISTREELs du variogramme/covariogramme
  120.       SEGACT MEVOL1
  121.       KEVOL1=MEVOL1.IEVOLL(1)
  122.       SEGACT KEVOL1
  123.       MLREE1=KEVOL1.IPROGX
  124.       MLREE2=KEVOL1.IPROGY
  125.  
  126. C     Rangement des donnees dans un segment SKRIG
  127.       MM=M
  128.       NN=N
  129.       NC=MLMOT1.MOTS(/2)
  130.       SEGINI SKRIG
  131.       DO K=1,NC
  132.           MOT1=MLMOT1.MOTS(K)
  133.           DO I=1,MM
  134.               IP1=IPT1.NUM(1,I)
  135.               CALL EXTRA9(MCHPO1,IP1,MOT1,IFOUR,.FALSE.,XFLOT,IRET)
  136.               MCOO(I,K)=XFLOT
  137.           ENDDO
  138.           DO I=1,NN
  139.               IP2=IPT2.NUM(1,I)
  140.               CALL EXTRA9(MCHPO2,IP2,MOT1,IFOUR,.FALSE.,XFLOT,IRET)
  141.               CCOO(I,K)=XFLOT
  142.           ENDDO
  143.       ENDDO
  144.       DO I=1,MM
  145.           IP1=IPT1.NUM(1,I)
  146.           CALL EXTRA9(MCHPO1,IP1,COMP,IFOUR,.FALSE.,XFLOT,IRET)
  147.           MVAL(I,1)=XFLOT
  148.       ENDDO
  149.  
  150.  
  151.  
  152. C----------------------------------------------------------------------C
  153. C                         CALCUL DES DISTANCES                         C
  154. C----------------------------------------------------------------------C
  155.  
  156. C     Entre les points de mesure 2 a 2
  157.       DO I=1,MM
  158.           DO J=I,MM
  159. C             Distance euclidienne
  160.               DIJ=XZERO
  161.               DO K=1,NC
  162.                   XI=MCOO(I,K)
  163.                   XJ=MCOO(J,K)
  164.                   DIJ=DIJ+((XI-XJ)**2)
  165.               ENDDO
  166.               DIJ=SQRT(DIJ)
  167.               DM(I,J)=DIJ
  168.               DM(J,I)=DIJ
  169.           ENDDO
  170.       ENDDO
  171.  
  172. C     Entre les points cibles et les points de mesure
  173.       DO I=1,MM
  174.           DO J=1,NN
  175. C             Distance euclidienne
  176.               DIJ=XZERO
  177.               DO K=1,NC
  178.                   XI=MCOO(I,K)
  179.                   XJ=CCOO(J,K)
  180.                   DIJ=DIJ+((XI-XJ)**2)
  181.               ENDDO
  182.               DIJ=SQRT(DIJ)
  183.               DC(I,J)=DIJ
  184.           ENDDO
  185.       ENDDO
  186.  
  187.  
  188.  
  189. C----------------------------------------------------------------------C
  190. C                      CALCUL DES (CO)VARIANCES :                      C
  191. C                  MATRICE ET VECTEURS SECOND MEMBRES                  C
  192. C----------------------------------------------------------------------C
  193.  
  194. C     Entre les points de mesure 2 a 2
  195.       DO I=1,MM
  196.           DO J=I,MM
  197.               DIJ=DM(I,J)
  198. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  199.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  200.               IF (IRET.NE.1) GOTO 999
  201.               VM(I,J)=VIJ
  202.               VM(J,I)=VIJ
  203.           ENDDO
  204.       ENDDO
  205. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  206.       DO J=1,MM
  207.           VM(MM+1,J)=XUN
  208.           VM(J,MM+1)=XUN
  209.       ENDDO
  210.       VM(MM+1,MM+1)=XZERO
  211.  
  212. C     Entre les points cibles et les points de mesure
  213.       DO I=1,MM
  214.           DO J=1,NN
  215.               DIJ=DC(I,J)
  216. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  217.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  218.               IF (IRET.NE.1) GOTO 999
  219.               VC(I,J)=VIJ
  220.               WC(I,J)=VIJ
  221.           ENDDO
  222.       ENDDO
  223. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  224.       DO J=1,NN
  225.           VC(MM+1,J)=XUN
  226.           WC(MM+1,J)=XUN
  227.       ENDDO
  228.  
  229.  
  230.  
  231. C----------------------------------------------------------------------C
  232. C                    RESOLUTION DU SYSTEME LINEAIRE                    C
  233. C----------------------------------------------------------------------C
  234.  
  235. C     Resolution du systeme VM * WC = VC a l'aide du solveur de LAPACK
  236.       CALL DGESV(MM+1,NN,VM,MM+1,IPIV,WC,MM+1,INFO)
  237. C     Sortie en cas d'erreur
  238.       IF (INFO.NE.0) THEN
  239.           PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO
  240.           PRINT*,'CHECK THE LINEAR SYSTEM'
  241.           CALL ERREUR(26)
  242.           RETURN
  243.       ENDIF
  244.  
  245.  
  246.  
  247. C----------------------------------------------------------------------C
  248. C         INTERPOLATION LINEAIRE DES VALEURS AUX POINTS CIBLES         C
  249. C----------------------------------------------------------------------C
  250.  
  251. C     Estimation de la fonction
  252.       JG=NN
  253.       SEGINI MLREE3
  254.       DO J=1,NN
  255.           XVAL=XZERO
  256.           DO I=1,MM
  257.               XVAL=XVAL+(WC(I,J)*MVAL(I,1))
  258.           ENDDO
  259.           MLREE3.PROG(J)=XVAL
  260.       ENDDO
  261.  
  262. C     Variance d'estimation
  263.       JG=NN
  264.       SEGINI MLREE4
  265.       DO J=1,NN
  266. C         Si on a utilise les variances, produit scalaire simple
  267.           IF (IVAR.EQ.1) THEN
  268.               XVAL=XZERO
  269.               DO I=1,MM+1
  270.                     XVAL=XVAL+(WC(I,J)*VC(I,J))
  271.               ENDDO
  272. C         Si on a utilise les covariances, petite soustraction a faire
  273.           ELSEIF (IVAR.EQ.2) THEN
  274.               CALL INTER5(XZERO,MLREE1,MLREE2,XSIG,0,0,1,IRET)
  275.               IF (IRET.NE.1) GOTO 999
  276.               XVAL=XSIG
  277.               DO I=1,MM+1
  278.                     XVAL=XVAL-(WC(I,J)*VC(I,J))
  279.               ENDDO
  280.           ENDIF
  281.           MLREE4.PROG(J)=XVAL
  282.       ENDDO
  283.  
  284.  
  285.  
  286. C----------------------------------------------------------------------C
  287. C                        ECRITURE DES RESULTATS                        C
  288. C----------------------------------------------------------------------C
  289.  
  290. C     Resultats sous forme de CHPOINT
  291. C     Variance d'estimation
  292.       CALL ECROBJ('LISTREEL',MLREE4)
  293.       CALL ECRCHA(COMP)
  294.       CALL ECROBJ('MAILLAGE',IPT2)
  295.       CALL MANUCH
  296.       CALL LIROBJ('CHPOINT',MCHPO4,1,IRETOU)
  297.       SEGACT MCHPO4*MOD
  298.       MCHPO4.JATTRI(1)=MCHPO2.JATTRI(1)
  299.       CALL ECROBJ('CHPOINT',MCHPO4)
  300. C     Estimation de la fonction
  301.       CALL ECROBJ('LISTREEL',MLREE3)
  302.       CALL ECRCHA(COMP)
  303.       CALL ECROBJ('MAILLAGE',IPT2)
  304.       CALL MANUCH
  305.       CALL LIROBJ('CHPOINT',MCHPO3,1,IRETOU)
  306.       SEGACT MCHPO3*MOD
  307.       MCHPO3.JATTRI(1)=MCHPO2.JATTRI(1)
  308.       CALL ECROBJ('CHPOINT',MCHPO3)
  309.  
  310. C     On fait le menage avant de partir
  311.       SEGDES MLMOT1
  312.       SEGDES MEVOL1
  313.       SEGDES KEVOL1
  314.       SEGSUP MLREE3,MLREE4
  315.       SEGDES MCHPO1
  316.       SEGDES MCHPO2,MCHPO3,MCHPO4
  317.  
  318. C     Et c'est fini !
  319.       RETURN
  320.  
  321. C     En cas d'erreur
  322. 999   CONTINUE
  323.       CALL ERREUR(26)
  324.       RETURN
  325.  
  326.       END
  327.  
  328.  
  329.  

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