Télécharger la plaquette Cast3M


fr en

Télécharger moca.eso

Retour à la liste

Numérotation des lignes : 

moca

  1. C MOCA      SOURCE    SP204843  26/02/16    21:15:07     12477          
  2.       subroutine moca
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8(A-H,O-Z)
  5.  
  6. -INC PPARAM
  7. -INC CCOPTIO
  8. -INC CCREEL 
  9. -INC SMLREEL
  10. -INC SMLOBJE
  11.       pointeur mlree4.mlreel, mlree5.mlreel
  12.       character*4 mot(1)
  13.       LOGICAL     LOK1
  14.       segment itra
  15.         integer ipoi(m),npara,ibon(m)
  16.       endsegment
  17.       segment itrav
  18.         real*8  a(n,n),b(n),c(n),bb(n,n)
  19.         integer is(n)
  20.       endsegment
  21.       data MOT /'POID'/
  22. * lecture de la valeur de depart des paramatres
  23.       call lirobj ('LISTREEL',mlpar,1,iretou)
  24.       mlreel=mlpar
  25.       segact mlreel
  26.       nparr = prog(/1)
  27. * lecture de la valeur theorique
  28.       call lirobj('LISTREEL',mlreel,1,iretou)
  29.       if(ierr.ne.0) return
  30.       segact mlreel
  31.       npmes=prog(/1)
  32. * lecture de la valeur de depart
  33.       call lirobj('LISTREEL',mlree1,1,iretou)
  34.       if(ierr.ne.0) return
  35.       segact mlree1
  36.       if(mlree1.prog(/1).ne.npmes) then
  37.       write(ioimp,*)' erreur 0 npmes mlree1.prog', npmes,mlree1.prog(/1)
  38.        call erreur (19)
  39.        return
  40.       endif
  41. * lecture des derivees par rapport au n parametres
  42. *     Cas de la donnee d'un LISTOBJE
  43.       CALL LIROBJ('LISTOBJE',ILOBJ1,0,IRET)
  44.       IF (IRET.NE.0) THEN
  45.         MLOBJE = ILOBJ1
  46.         SEGACT,MLOBJE
  47.         IF (TYPOBJ.NE.'LISTREEL') THEN
  48.           MOTERR(1:8)='LISTREEL'
  49.           CALL ERREUR(1156)
  50.           RETURN
  51.         ENDIF
  52.         NOBJ1 = LISOBJ(/1)
  53.         M = NOBJ1 
  54.         SEGINI,ITRA
  55.         NPARA = 0
  56.         IF (NOBJ1.NE.NPARR) THEN
  57.           INTERR = NPARR
  58.           CALL ERREUR(1157)
  59.           RETURN
  60.         ENDIF
  61.         ZPREC1 = 10.*XZPREC
  62.         DO IO1=1,NOBJ1
  63.           MLREE2 = LISOBJ(IO1)
  64.           SEGACT,MLREE2
  65.           IF (MLREE2.PROG(/1).NE.NPMES) THEN
  66.             INTERR = NPMES
  67.             CALL ERREUR(1158)
  68.             RETURN
  69.           ENDIF
  70.           LOK1 = .FALSE.
  71.           DO IL1=1,MLREE2.PROG(/1)
  72.             XVAL1 = ABS(MLREE2.PROG(IL1))
  73.             LOK1 = XVAL1.GT.((1.-ZPREC1)*XVAL1)
  74.             IF (LOK1) GOTO 30
  75.           ENDDO
  76.  30       CONTINUE
  77.           IF (LOK1) THEN
  78.             IBON(IO1) = 1
  79.             NPARA = NPARA+1
  80.             IPOI(NPARA) = MLREE2
  81.           ENDIF
  82.         ENDDO
  83.       ELSE 
  84. *     Cas de la donnee de plusieurs LISTREEL
  85.         m=100
  86.         segini itra
  87.         NPARA = 0
  88. *        write(ioimp,*) ' nb de parametres ',nparr
  89.         do 1 i=1,nparr
  90.          call lirobj('LISTREEL',mlree2,1,iretou)
  91.          if(ierr.ne.0) then
  92.           write(ioimp,*)  ' erreur 1'
  93.            call erreur (19)
  94.            return
  95.          endif
  96.          segact mlree2
  97.          do iy=1,mlree2.prog(/1)
  98.           if( mlree2.prog(iy).ne.0.d0) go to 10
  99.          enddo
  100.          go to 1
  101.    10    continue
  102.          ibon(i)=1
  103.          if(mlree2.prog(/1).ne.npmes) then
  104.       write(ioimp,*)' erreur 2  npmes mlree2.prog',npmes,mlree2.prog(/1)
  105.            call erreur (19)
  106.            return
  107.          endif
  108.          npara=npara+1
  109.          ipoi(npara)=mlree2
  110.    1    continue
  111.       ENDIF 
  112.       if (npara.eq.0) then
  113.         write(ioimp,*)  ' erreur 3  npara ',npara
  114.         call erreur (21)
  115.         return
  116.       endif
  117. * lecture de l'option POIDS
  118.       call lirmot(mot,1,itrou,0)
  119.       if( itrou.eq.0) then
  120.         jg=npmes
  121.          isup=1
  122.         segini mlree5
  123.         do io=1,npmes
  124.           mlree5.prog(io)=1.d0
  125.         enddo
  126.       else
  127. * lecture du poids
  128.         call lirobj('LISTREEL',mlree5,1,iretou)
  129.         if(ierr.ne.0)return
  130.         segact mlree5
  131.         isup=0
  132.       endif
  133.       n = npara
  134.       segini itrav
  135.       do 3 i=1,npmes
  136.       do 4 j=1,npara
  137.        mlree3=ipoi(j)
  138.        x1=mlree3.prog(i)*2. *mlree5.prog(i)*mlree5.prog(i)
  139.        do 5 k=1,npara
  140.         mlree4=ipoi(k)
  141.         xk=mlree4.prog(i)
  142.         a(j,k)=a(j,k)+x1*xk
  143.    5   continue
  144.        b(j) = b(j) + x1 * ( mlreel.prog(i)-mlree1.prog(i))
  145.    4  continue
  146.    3  continue
  147. *
  148. **  appel a l'inversion
  149. *
  150. *      write(6,*) ' avant inversion'
  151. *      write(6,*) ( a(1,it),it=1,npara)
  152. *      write(6,*) ( a(2,it),it=1,npara)
  153. *      write(6,*) ( a(3,it),it=1,npara)
  154.  
  155.       diamax=0.
  156.       do io=1,npara
  157.         if( abs(a(io,io)).gt.diamax) diamax=abs(a(io,io))
  158.       enddo
  159.       eps = diamax / 1.e10
  160.       call inver (A,npara,icrit,c,is,eps)
  161.       if(icrit.ne.0) then
  162.        write(ioimp,*) '***** Inversion systeme impossible'
  163.        call erreur(21)
  164.        return
  165.       endif
  166. *      write(6,*) ' apres inversion'
  167. *      write(6,*) ( a(1,it),it=1,npara)
  168. *      write(6,*) ( a(2,it),it=1,npara)
  169. *      write(6,*) ( a(3,it),it=1,npara)
  170.  
  171. *      write(6,*) ' b '
  172. *      write(6,*) ( b(it),it=1,npara)
  173.       if(isup.eq.1) then
  174.        segsup mlree5
  175.       else
  176.        segdes mlree5
  177.       endif
  178.       jg=npara
  179.       segini mlree5
  180.       do 50 io=1,npara
  181.       mlree5.prog(io)=0.d0
  182.       do 6 iu=1,npara
  183.       mlree5.prog(io)=mlree5.prog(io)+a(io,iu)*b(iu)
  184.    6  continue
  185.   50  continue
  186. *      write(6,*) 'valeur trouvées dans moca'
  187. *      write(6,*) (mlree5.prog(io),io=1,npara)
  188.       segdes mlreel,mlree1
  189.       do io=1,npara
  190.       mlree3=ipoi(io)
  191.       segdes mlree3
  192.       enddo
  193.       jg=nparr
  194.       segini mlreel
  195.       ia=1
  196.       mlree2=mlpar
  197.       do ib=1,nparr
  198.       if(ibon(ib).eq.1)then
  199.       prog(ib)=mlree5.prog(ia) + mlree2.prog(ib)
  200.       ia=ia+1
  201.        else
  202.        prog(ib)=mlree2.prog(ib)
  203.       endif
  204.       enddo
  205.       segdes mlreel,mlree2
  206.       segsup mlree5
  207.       call ecrobj('LISTREEL',mlreel)
  208.       segsup itrav,itra
  209.       return
  210.       end
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  

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