Esope source files

   1 : C MTCLSD    SOURCE    CHAT      05/01/13    01:53:48     5004           
   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,sig(3)
  14 : 
  14 :       call zzero(amat ,16)
  14 :       call zzero(aamat,16)
  14 :       call zzero(eemat,16)
  14 :       call zzero(gmat ,16)
  14 :       call zzero(kmat ,16)
  14 :       call zzero(vecm,4)
  14 :       call zzero(vecn,4)
  14 :       call zzero(vecnaux,4)
  14 :       call zzero(vecn,4)
  14 :       call zzero(Avecm,4)
  14 :       call zzero(vecntA,4)
  14 :       call zzero(sig,3)
  15 :       void=0.D0
  16 :       ndimv=1
  17 :       ndims=3
  17 :       call der_enerelas_dpral(x,sig,nmodel)
  18 : c m,n,A
  19 :       if (ndimx.eq.3) then
  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 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
  35 :       call zzero(amat,16 )
  36 : c A=I+l*AA
  37 :       do i=1,ndimx
  38 :        amat(i+(i-1)*ndimx)=1.D0
  39 :        do j=1,ndimx
  40 :          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  41 :      .                       lam*aamat(i+(j-1)*ndimx)
  42 :        enddo
  43 :       enddo
  44 : c G = A-1
  44 :       call DescLU(Amat,ndimx)
  44 :       call LUinv(Amat,Gmat,ndimx)
  45 : c g = nt*A-1*m
  46 :       g=0.D0
  47 :       do i=1,ndimx
  48 :        do j=1,ndimx
  49 :           g=g+vecn(i)*Gmat(i+ndimx*(j-1))*vecm(j)
  50 :        enddo
  51 :       enddo
  52 : c G*m; nt*G
  53 :       do i=1,ndimx
  54 :        do j=1,ndimx
  55 :          Avecm(i)=Avecm(i)+Gmat(i+ndimx*(j-1))*vecm(j)
  56 :          vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndimx*(i-1))
  57 :        enddo
  58 :       enddo
  59 : c K = G-(G*m * nt*G)/g
  60 :       do i=1,ndimx
  61 :        do j=1,ndimx
  62 :           Kmat(i,j)=Gmat((j-1)*ndimx+i)-Avecm(i)*vecntA(j)/g
  63 :        enddo
  64 :       enddo
  65 : c EK= E*K
  65 :       call zzero(ekmat,ndimk**2)
  66 :       do i=1,3
  67 :        do j=1,3
  68 :          ekmat(i,j)=0.D0
  69 :          do k=1,ndimx
  70 :           ekmat(i,j)=ekmat(i,j)+
  71 :      .               eemat(i+(k-1)*ndimx)*kmat(k,j)
  72 :          enddo 
  73 :        enddo
  74 :       enddo
  75 : 
  76 :       if (nescri.eq.1)
  77 :      .   write(ues,999)' EKmat_dp ',((ekmat(i,j),j=1,3),i=1,3)
  78 : 999   format(2x,A12,/,3(3(1x,E12.6),/))
  79 :       return
  80 :       end
  81 :  
  82 :  
  83 :