Télécharger la plaquette Cast3M


fr en

Télécharger mtclsd.eso

Retour à la liste

Numérotation des lignes : 

mtclsd

  1. C MTCLSD    SOURCE    PV        22/04/19    16:18:08     11344          
  2. CCC
  3. C **********************************************************************
  4. CCC
  5.       SUBROUTINE MTC_LS_DPRAL(EKMAT,NDIMK,X,NDIMX,LAM,
  6.      .                        nmodel,nescri,ues,nnumer,deltax)
  7.  
  8.       IMPLICIT INTEGER(I-N)
  9.       integer ndimk,ndimx,nmodel,nescri,ues,nnumer
  10.       integer i,j,k,ndims,ndimv
  11.       real*8  EKmat(ndimk,ndimk),x(ndimx),lam,deltax
  12.       real*8  amat(16),aamat(16),eemat(16),vecm(4),vecnaux(4),vecn(4),
  13.      .        gmat(16),Avecm(4),vecntA(4),kmat(4,4),g,void(1),sig(3)
  14.  
  15.       call zzero(amat ,16)
  16.       call zzero(aamat,16)
  17.       call zzero(eemat,16)
  18.       call zzero(gmat ,16)
  19.       call zzero(kmat ,16)
  20.       call zzero(vecm,4)
  21.       call zzero(vecn,4)
  22.       call zzero(vecnaux,4)
  23.       call zzero(vecn,4)
  24.       call zzero(Avecm,4)
  25.       call zzero(vecntA,4)
  26.       call zzero(sig,3)
  27.       void(1)=0.D0
  28.       ndimv=1
  29.       ndims=3
  30.       call der_enerelas_dpral(x,sig,nmodel)
  31. c m,n,A
  32.       if (ndimx.eq.3) then
  33.        call vflsig(sig,ndims,void,ndimv,vecm,nmodel)
  34.        call vyisig(sig,ndims,void,ndimv,vecnaux,nmodel)
  35.        call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
  36.       else if (ndimx.eq.4) then
  37.        call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
  38.        call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
  39.        call vyisig(sig,ndims,x(ndimx),ndimv,vecnaux,nmodel)
  40.        call vyivar(sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel)
  41.        call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  42.       endif
  43. c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  44.       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  45. c n^T=n^T*E
  46. c AA=A*E
  47.       do i=1,ndimx
  48.        vecn(i)=0.D0
  49.        do j=1,ndimx
  50.          vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx)
  51.          aamat(i+(j-1)*ndimx)=0.D0
  52.          do k=1,ndimx
  53.            aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
  54.      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  55.          enddo
  56.        enddo
  57.       enddo
  58.       call zzero(amat,16 )
  59. c A=I+l*AA
  60.       do i=1,ndimx
  61.        amat(i+(i-1)*ndimx)=1.D0
  62.        do j=1,ndimx
  63.          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  64.      .                       lam*aamat(i+(j-1)*ndimx)
  65.        enddo
  66.       enddo
  67. c G = A-1
  68.       call DescLU(Amat,ndimx)
  69.       call LUinv(Amat,Gmat,ndimx)
  70. c g = nt*A-1*m
  71.       g=0.D0
  72.       do i=1,ndimx
  73.        do j=1,ndimx
  74.           g=g+vecn(i)*Gmat(i+ndimx*(j-1))*vecm(j)
  75.        enddo
  76.       enddo
  77. c G*m; nt*G
  78.       do i=1,ndimx
  79.        do j=1,ndimx
  80.          Avecm(i)=Avecm(i)+Gmat(i+ndimx*(j-1))*vecm(j)
  81.          vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndimx*(i-1))
  82.        enddo
  83.       enddo
  84. c K = G-(G*m * nt*G)/g
  85.       do i=1,ndimx
  86.        do j=1,ndimx
  87.           Kmat(i,j)=Gmat((j-1)*ndimx+i)-Avecm(i)*vecntA(j)/g
  88.        enddo
  89.       enddo
  90. c EK= E*K
  91.       call zzero(ekmat,ndimk**2)
  92.       do i=1,3
  93.        do j=1,3
  94.          ekmat(i,j)=0.D0
  95.          do k=1,ndimx
  96.           ekmat(i,j)=ekmat(i,j)+
  97.      .               eemat(i+(k-1)*ndimx)*kmat(k,j)
  98.          enddo
  99.        enddo
  100.       enddo
  101.  
  102.       if (nescri.eq.1)
  103.      .   write(ues,999)' EKmat_dp ',((ekmat(i,j),j=1,3),i=1,3)
  104. 999   format(2x,A12,/,3(3(1x,E12.6),/))
  105.       return
  106.       end
  107.  
  108.  
  109.  
  110.  
  111.  

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