1 : C MTCLDD SOURCE PV 07/11/23 21:17:54 5978 2 : CCC 3 : C ********************************************************************** 4 : CCC 5 : SUBROUTINE MTC_LS_DPRAL_DENSI(EKMAT,NDIMK,X,NDIMX,LAM, 6 : . xdensi,xmat, 7 : . nmodel,nescri,ues,nnumer,deltax) 8 : 9 : IMPLICIT INTEGER(I-N) 10 : IMPLICIT REAL*8 (a-h,o-z) 11 : integer ndimk,ndimx,nmodel,nescri,ues,nnumer 12 : integer i,j,k,ndims,ndimv 13 : real*8 EKmat(ndimk,ndimk),x(ndimx),lam,deltax 14 : real*8 amat(16),aamat(16),eemat(16),vecm(4),vecnaux(4),vecn(4), 15 : . gmat(16),Avecm(4),vecntA(4),kmat(4,4),g,void,sig(3) 16 : real*8 xdensi,dxdensi,xdensi2,xmat(*),resu,resuaux 17 : real*8 Mmat(4,3),kmat2(4,4),kmat3(4,4),ekmat2(3,3),ekmat3(3,3), 18 : . vecmaux(4) 19 : 19 : call zzero(amat ,16) 19 : call zzero(aamat,16) 19 : call zzero(eemat,16) 19 : call zzero(gmat ,16) 19 : call zzero(kmat ,16) 19 : call zzero(vecm,4) 19 : call zzero(vecn,4) 19 : call zzero(vecnaux,4) 19 : call zzero(vecn,4) 19 : call zzero(Avecm,4) 19 : call zzero(vecntA,4) 19 : call zzero(sig,3) 20 : 20 : call zzero(mmat,12) 20 : call zzero(kmat2,16) 20 : call zzero(kmat3,16) 20 : call zzero(ekmat2,9) 20 : call zzero(ekmat3,9) 20 : call zzero(vecmaux,4) 21 : 22 : void=0.D0 23 : ndimv=1 24 : ndims=3 24 : call der_enerelas_dpral(x,sig,nmodel) 25 : c m,n,A 26 : if (ndimx.eq.3) then 26 : call yieldd(sig,ndims,void,ndimv,resu,nmodel) 26 : call vflsig(sig,ndims,void,ndimv,vecm,nmodel) 26 : call vyisig(sig,ndims,void,ndimv,vecnaux,nmodel) 26 : call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel) 27 : else if (ndimx.eq.4) then 27 : call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel) 27 : call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel) 27 : call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel) 27 : call vyisig(sig,ndims,x(ndimx),ndimv,vecnaux,nmodel) 27 : call vyivar(sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel) 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) 29 : call der2_enerelas_dpral(x,eemat,ndimx,nmodel) 30 : c n^T=n^T*E 31 : c AA=A*E 32 : do i=1,ndimx 33 : vecn(i)=0.D0 34 : do j=1,ndimx 35 : vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx) 36 : aamat(i+(j-1)*ndimx)=0.D0 37 : do k=1,ndimx 38 : aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+ 39 : . amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx) 40 : enddo 41 : enddo 42 : enddo 42 : call zzero(amat,16 ) 43 : c A=I+l*AA 44 : do i=1,ndimx 45 : amat(i+(i-1)*ndimx)=1.D0 46 : do j=1,ndimx 47 : amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+ 48 : . lam*aamat(i+(j-1)*ndimx) 49 : enddo 50 : enddo 51 : c G = A-1 51 : call DescLU(Amat,ndimx) 51 : call LUinv(Amat,Gmat,ndimx) 52 : c g = nt*A-1*m 53 : g=0.D0 54 : do i=1,ndimx 55 : do j=1,ndimx 56 : g=g+vecn(i)*Gmat(i+ndimx*(j-1))*vecm(j) 57 : enddo 58 : enddo 59 : c G*m; nt*G 60 : do i=1,ndimx 61 : do j=1,ndimx 62 : Avecm(i)=Avecm(i)+Gmat(i+ndimx*(j-1))*vecm(j) 63 : vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndimx*(i-1)) 64 : enddo 65 : enddo 66 : c K = G-(G*m * nt*G)/g 67 : do i=1,ndimx 68 : do j=1,ndimx 69 : Kmat(i,j)=Gmat((j-1)*ndimx+i)-Avecm(i)*vecntA(j)/g 70 : enddo 71 : enddo 72 : CCCC 73 : C densi 74 : CCCC 75 : c Vecm y f perturbados 76 : 77 : c write(399,*) 78 : c dxrel = 0.1 79 : c do i = 1,10 80 : c dx = xdensi * (dxrel**i) 81 : c xdensi2 = xdensi + dx 81 : c call carac_mate_rhmc_densi(XMAT,xdensi2) 81 : c call yieldd(sig,ndims,void,ndimv,resuaux,nmodel) 81 : c call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel) 82 : c write(399,*)dxrel**i,(resuaux-resu)/dx 83 : c write(399,*)' ',((vecmaux(j)-vecm(j))/dx,j=1,3) 83 : c call carac_mate_rhmc_densi(XMAT,xdensi) 84 : c enddo 85 : c do i = 1,10 86 : c dx = - xdensi * (dxrel**i) 87 : c xdensi2 = xdensi + dx 87 : c call carac_mate_rhmc_densi(XMAT,xdensi2) 87 : c call yieldd(sig,ndims,void,ndimv,resuaux,nmodel) 87 : c call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel) 88 : c write(399,*)-dxrel**i,(resuaux-resu)/dx 89 : c write(399,*)' ',((vecmaux(j)-vecm(j))/dx,j=1,3) 89 : c call carac_mate_rhmc_densi(XMAT,xdensi) 90 : c enddo 91 : CC 92 : dxrel = 0.1D0 93 : dx = (1.D0 + xdensi) * (dxrel**6) 94 : xdensi2 = xdensi + dx 94 : c call carac_mate_rhmc_densi(XMAT,xdensi2) 94 : call carac_mate_densi(XMAT,xdensi2,nmodel) 95 : if (ndimx.eq.3) then 95 : call yieldd(sig,ndims,void,ndimv,resuaux,nmodel) 95 : call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel) 96 : else if (ndimx.eq.4) then 96 : call yieldd(sig,ndims,x(ndimx),ndimv,resuaux,nmodel) 96 : call vflsig(sig,ndims,x(ndimx),ndimv,vecmaux,nmodel) 96 : call vflvar(sig,ndims,x(ndimx),ndimv,vecmaux(ndimx),nmodel) 97 : endif 98 : xdensi2 = xdensi - dx 98 : c call carac_mate_rhmc_densi(XMAT,xdensi2) 98 : call carac_mate_densi(XMAT,xdensi2,nmodel) 99 : if (ndimx.eq.3) then 99 : call yieldd(sig,ndims,void,ndimv,resu,nmodel) 99 : call vflsig(sig,ndims,void,ndimv,vecm,nmodel) 100 : else if (ndimx.eq.4) then 100 : call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel) 100 : call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel) 100 : call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel) 101 : endif 102 : dx = dx * 2.D0 102 : c call carac_mate_rhmc_densi(XMAT,xdensi) 102 : call carac_mate_densi(XMAT,xdensi,nmodel) 103 : c Mmat = derivada de vecm (3 veces) 104 : do i=1,ndimx 105 : Mmat(i,1)=(vecmaux(i)-vecm(i))/dx 106 : Mmat(i,2)=Mmat(i,1) 107 : Mmat(i,3)=Mmat(i,1) 108 : enddo 109 : c Kmat2 = Kmat * Mmat 110 : c Kmat3 = Avecm/g (3 veces) 111 : do i=1,ndimx 112 : do j=1,3 113 : do k=1,ndimx 114 : Kmat2(i,j)=Kmat2(i,j)+Kmat(i,k)*Mmat(k,j) 115 : enddo 116 : Kmat3(i,j)=Avecm(i)/g 117 : enddo 118 : enddo 119 : c auxiliares 120 : aux1= xdensi*lam 121 : aux2= xdensi*(resuaux-resu)/dx 122 : c EKmat = E*Kmat ; EKmat2 = E*Kmat2 ; EKmat3 = E*Kmat3 123 : c EKmat = EKmat + aux1*EKmat2 + aux2*EKmat3 123 : call zzero(ekmat,ndimk**2) 124 : do i=1,3 125 : do j=1,3 126 : do k=1,ndimx 127 : ekmat(i,j) =ekmat(i,j)+ 128 : . eemat(i+(k-1)*ndimx)*kmat(k,j) 129 : ekmat2(i,j)=ekmat2(i,j)+ 130 : . eemat(i+(k-1)*ndimx)*kmat2(k,j) 131 : ekmat3(i,j)=ekmat3(i,j)+ 132 : . eemat(i+(k-1)*ndimx)*kmat3(k,j) 133 : enddo 134 : ekmat(i,j)=ekmat(i,j)+aux1*ekmat2(i,j)+aux2*ekmat3(i,j) 135 : enddo 136 : enddo 137 : 138 : if (nescri.eq.1) 139 : . write(ues,999)' EKmat_dp ',((ekmat(i,j),j=1,3),i=1,3) 140 : 999 format(2x,A12,/,3(3(1x,E12.6),/)) 141 : return 142 : end 143 : 144 : 145 : 146 :
© Cast3M 2003 - All rights reserved.
Disclaimer