Télécharger la plaquette Cast3M


fr en

Télécharger di2lev.eso

Retour à la liste

Numérotation des lignes : 

di2lev

  1. C DI2LEV    SOURCE    PV        22/04/19    16:18:01     11344          
  2.  
  3.       SUBROUTINE DIRECTION_2LEVELS(X,R,DX,N,NDIMX,NMODEL,SIG,VECM,
  4.      .                             nnumer,deltax,lam)
  5.       IMPLICIT INTEGER(I-N)
  6.       integer n,i,j,k,ndims,ndimv,nmodel,nnumer,ndimx
  7.       real*8  r(n),x(n),dx(n),deltax,sig(3),vecm(4)
  8.       real*8  vecnaux(4),vecn(4),amat(16),vtLUiw,void(1)
  9.       real*8  aamat(16),eemat(16),lam
  10.       integer             augla
  11.       real*8                    c
  12.       common /auglagrang1/ augla
  13.       common /auglagrang2/ c
  14.       real*8  bbmat(16),res
  15.       call zzero(amat ,16)
  16.       call zzero(aamat,16)
  17.       call zzero(bbmat,16)
  18.       call zzero(eemat,16)
  19.       void(1)=0.D0
  20.       ndimv=1
  21.       ndims=3
  22. c n
  23. c A
  24.       if (ndimx.eq.3) then
  25.        call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
  26.       else if (ndimx.eq.4) then
  27.        call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  28.       endif
  29. c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  30.       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  31. c AA=A*E
  32.       do i=1,ndimx
  33.        do j=1,ndimx
  34.          aamat(i+(j-1)*ndimx)=0.D0
  35.          do k=1,ndimx
  36.           aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
  37.      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  38.          enddo
  39.        enddo
  40.       enddo
  41. **
  42.       if ((augla.eq.1)) then
  43. c residuo y v flujo
  44.       if (ndimx.eq.3) then
  45.        call yieldd(sig,ndims,void,ndimv,res,nmodel)
  46.        call vyisig(sig,ndims,void,ndimv,vecn,nmodel)
  47.       else if (ndimx.eq.4) then
  48.        call yieldd(sig,ndims,x(ndimx),ndimv,res,nmodel)
  49.        call vyisig(sig,ndims,x(ndimx),ndimv,vecn,nmodel)
  50.        call vyivar(sig,ndims,x(ndimx),ndimv,vecn(ndimx),nmodel)
  51.       endif
  52. c BB=n*n^T*E
  53.       do i=1,ndimx
  54.        do j=1,ndimx
  55.          bbmat(i+(j-1)*ndimx)=0.D0
  56.            do k=1,ndimx
  57.             bbmat(i+(j-1)*ndimx)=bbmat(i+(j-1)*ndimx)+
  58.      .                         vecn(i)*vecn(k)*eemat(k+(j-1)*ndimx)
  59.            enddo
  60.        enddo
  61.       enddo
  62.       endif
  63. **
  64.       call zzero(amat,16 )
  65. c A=I+l*AA
  66.       do i=1,ndimx
  67.        amat(i+(i-1)*ndimx)=1.D0
  68.        do j=1,ndimx
  69. c         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  70. c     .                       lam*aamat(i+(j-1)*ndimx)
  71. **
  72.          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  73.      .                       ABS(lam+c*res)*aamat(i+(j-1)*ndimx)+
  74.      .                       c*bbmat(i+(j-1)*ndimx)
  75. **
  76.        enddo
  77.       enddo
  78.       call DescLU(Amat,ndimx)
  79.       call LUiw(Amat,r,dx,ndimx)
  80.       do i=1,ndimx
  81.         dx(i)=-1.D0*dx(i)
  82.       enddo
  83.       return
  84.       end
  85.  
  86.  
  87.  
  88.  
  89.  
  90.  
  91.  

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