Esope source files

   1 : C ACT_LA    SOURCE    CHAT      05/01/12    21:16:48     5004           
   2 : 
   3 :       SUBROUTINE ACTUALIZA_LAMBDA(X,NDIMX,NMODEL,LAM,DLA,DLB)
   4 :       IMPLICIT INTEGER(I-N)
   5 :       integer n,i,j,k,ndims,ndimv,nmodel,nnumer,ndimx
   6 :       real*8  x(ndimx),deltax,sig(3),vecm(4)
   7 :       real*8  vecnaux(4),vecn(4),amat(16),vtLUiw,void
   8 :       real*8  aamat(16),eemat(16),dla,dlb,lam
   9 :       integer             augla
  10 :       real*8                    c
  11 :       common /auglagrang1/ augla
  12 :       common /auglagrang2/ c
  13 :       real*8  bbmat(16),res
  13 :       call zzero(amat ,16)
  13 :       call zzero(aamat,16)
  13 :       call zzero(bbmat,16)
  13 :       call zzero(eemat,16)
  14 :       void=0.D0
  15 :       ndimv=1
  16 :       ndims=3
  17 : c conseguir sig
  17 :       call der_enerelas_dpral(x,sig,nmodel)
  18 : c n, m , A
  19 :       if (ndimx.eq.3) then
  19 :        call yieldd(sig,ndims,void,ndimv,dla,nmodel)
  19 :        call vflsig(sig,ndims,void,ndimv,vecm,nmodel)
  19 :        call vyisig(sig,ndims,void,ndimv,vecnaux,nmodel)
  19 :        call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
  20 :       else if (ndimx.eq.4) then 
  20 :        call yieldd(sig,ndims,x(ndimx),ndimv,dla,nmodel)
  20 :        call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
  20 :        call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
  20 :        call vyisig(sig,ndims,x(ndimx),ndimv,vecnaux,nmodel)
  20 :        call vyivar(sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel)
  20 :        call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  21 :       endif
  22 : c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  22 :       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  23 : c n^T=n^T*E
  24 : c AA=A*E
  25 :       do i=1,ndimx
  26 :        vecn(i)=0.D0
  27 :        do j=1,ndimx
  28 :          vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx)
  29 :          aamat(i+(j-1)*ndimx)=0.D0
  30 :          do k=1,ndimx
  31 :           aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
  32 :      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  33 :          enddo 
  34 :        enddo
  35 :       enddo
  36 : **
  37 :       if ((augla.eq.1)) then 
  38 : c BB=n*n^T*E
  39 :       do i=1,ndimx
  40 :        do j=1,ndimx
  41 :          bbmat(i+(j-1)*ndimx)=0.D0
  42 :          do k=1,ndimx
  43 :           bbmat(i+(j-1)*ndimx)=bbmat(i+(j-1)*ndimx)+
  44 :      .                      vecnaux(i)*vecnaux(k)*eemat(k+(j-1)*ndimx)
  45 :          enddo 
  46 :        enddo
  47 :       enddo
  48 :       endif
  49 : **
  49 :       call zzero(amat,16 )
  50 : c A=I+l*AA
  51 :       do i=1,ndimx
  52 :        amat(i+(i-1)*ndimx)=1.D0
  53 :        do j=1,ndimx
  54 : c         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  55 : c     .                       lam*aamat(i+(j-1)*ndimx)
  56 : **
  57 :          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  58 :      .                       dabs(lam+c*dla)*aamat(i+(j-1)*ndimx)+
  59 :      .                       c*bbmat(i+(j-1)*ndimx)
  60 : **
  61 :        enddo
  62 :       enddo
  62 :       call DescLU(Amat,ndimx)
  63 :       dlb=-vtLUiw(vecn,Amat,vecm,ndimx)
  64 :       return
  65 :       end
  66 :  
  67 :  
  68 :