[cig-commits] r6262 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu Mar 15 12:38:25 PDT 2007
Author: willic3
Date: 2007-03-15 12:38:25 -0700 (Thu, 15 Mar 2007)
New Revision: 6262
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
Log:
Major fixes to Newton-Raphson solution for power-law.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f 2007-03-15 17:15:26 UTC (rev 6261)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/lnsrch.f 2007-03-15 19:38:25 UTC (rev 6262)
@@ -1,4 +1,4 @@
- SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func,rpar,
+ SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,fvec,stpmax,check,func,rpar,
& nrpar,ipar,nipar,funcv,ierr,errstrng)
c
c... routine to perform a line search along the direction given by p.
@@ -18,7 +18,7 @@
c
integer n,nrpar,nipar,ierr
integer ipar(nipar)
- double precision xold(n),fold,g(n),p(n),x(n),f,stpmax,func
+ double precision xold(n),fold,g(n),p(n),x(n),f,fvec(n),stpmax,func
double precision rpar(nrpar)
external func,funcv
logical check
@@ -35,7 +35,6 @@
INTEGER i
double precision a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2
double precision slope,sum,temp,test,tmplam,da
- double precision fvec(NP)
c
check=.false.
sum=dnrm2(n,p,ione)
@@ -47,7 +46,7 @@
test=zero
do i=1,n
temp=abs(p(i))/max(abs(xold(i)),one)
- if(temp.gt.test)test=temp
+ test=max(temp,test)
end do
alamin=TOLX/test
alam=one
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-15 17:15:26 UTC (rev 6261)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2007-03-15 19:38:25 UTC (rev 6262)
@@ -323,11 +323,13 @@
end
c
c
- subroutine jaccmp_6(scur,n,fvec,NP,fjac,rpar,nrpar,ipar,nipar,
- & ierr,errstrng)
+ subroutine jaccmp_6(scur,n,fvec,nfjsize,fjac,rpar,nrpar,ipar,
+ & nipar,ierr,errstrng)
c
c... subroutine to form the Jacobian used in finding the zero of the
c stress function.
+c** Note that in this version fjac is always assumed to be a symmetric
+c matrix stored in packed format.
c
include "implicit.inc"
c
@@ -340,11 +342,11 @@
c
c... subroutine arguments
c
- integer n,NP,nrpar,nipar,ierr
+ integer n,nfjsize,nrpar,nipar,ierr
c*** note that the dimension below is a bit kludgy and does not allow
c*** anything else to be put in the ipar array
integer ipar(nstr,nstr)
- double precision scur(n),fvec(n),fjac(NP,NP),rpar(nrpar)
+ double precision scur(n),fvec(n),fjac(nfjsize),rpar(nrpar)
character errstrng*(*)
c
c... local variables
@@ -352,7 +354,7 @@
integer i,j
ctest integer iddmat(nstr,nstr)
ctest equivalence(ipar(indiddmat),iddmat)
- double precision dtmp(nddmat),cmat(nddmat)
+ double precision cmat(nddmat)
double precision sdev(nstr),sinv1,steff
c
cdebug write(6,*) "Hello from jaccmp_6_f!"
@@ -361,22 +363,14 @@
c... get stress invariants and a copy of inverted elasticity matrix
c
call invar(sdev,sinv1,steff,scur)
- call dcopy(nddmat,rpar(inddmati),ione,dtmp,ione)
+ call dcopy(nddmat,rpar(inddmati),ione,fjac,ione)
c
c... if second deviatoric invariant is zero, use only elastic solution
c
if(steff.eq.zero) then
- call invsymp(nddmat,nstr,dtmp,ierr,errstrng)
- if(ierr.ne.izero) return
- do i=1,nstr
- do j=1,nstr
- fjac(i,j)=dtmp(ipar(i,j))
- end do
- end do
return
c
-c... otherwise, invert elastic matrix and augment it by the viscous
-c Jacobian.
+c... otherwise, augment inverted elastic matrix by the viscous Jacobian.
c
else
c
@@ -384,16 +378,7 @@
c
call dbds_6(cmat,nddmat,ipar,sdev,nstr,steff,
& rpar(indanpwr),rpar(indemhu),rpar(indalfap),rpar(inddeltp))
- call daxpy(nddmat,one,cmat,ione,dtmp,ione)
-c
-c... invert augmented matrix and store results in fjac
-c
- call invsymp(nddmat,nstr,dtmp,ierr,errstrng)
- do i=1,nstr
- do j=1,nstr
- fjac(i,j)=dtmp(ipar(i,j))
- end do
- end do
+ call daxpy(nddmat,one,cmat,ione,fjac,ione)
end if
return
end
@@ -548,7 +533,7 @@
c
integer i
double precision e,pr,anpwr,emhu,rmu
- double precision test0,rm,rmt,rmtdt
+ double precision test0,rmt,rmtdt
double precision sdev(nstr),betat(nstr),betatdt(nstr),sinv1,seff
double precision rpar(nrpar)
logical check
@@ -573,7 +558,6 @@
rmu=half*e/(one+pr)
tmax=big
rmt=deltp*(one-alfap)
- rm=sub*rmt
rmtdt=deltp*alfap
c
c... for first iteration, use stresses from previous step as initial
@@ -585,17 +569,18 @@
c
c... current total strain
call dcopy(nstr,ee,ione,rpar(indconst),ione)
+ call dscal(nstr,sub,rpar(indconst),ione)
c... inverted elasticity matrix times initial stresses
call elas_mat_6(rpar(inddmati),prop,iddmat,nprop,ierr,errstrng)
call invsymp(nddmat,nstr,rpar(inddmati),ierr,errstrng)
if(ierr.ne.izero) return
- if(test0.ne.zero) call dspmv("u",nstr,one,rpar(inddmati),state0,
+ if(test0.ne.zero) call dspmv("u",nstr,sub,rpar(inddmati),state0,
& ione,one,rpar(indconst),ione)
c... viscous strain from previous time step
- call daxpy(nstr,sub,state(13),ione,rpar(indconst),ione)
+ call daxpy(nstr,one,state(13),ione,rpar(indconst),ione)
c... contribution to viscous strain from stresses at beginning of step
call betacmp_6(state,betat,sdev,seff,sinv1,emhu,anpwr,nstr)
- call daxpy(nstr,rm,betat,ione,rpar(indconst),ione)
+ call daxpy(nstr,rmt,betat,ione,rpar(indconst),ione)
c
c... finish defining parameters for stress computation then call Newton
c routine to compute stresses.
@@ -673,7 +658,7 @@
c
c... viscous strain rate multiplier
c
- rmtdt=-rpar(indalfap)*rpar(inddeltp)
+ rmtdt=rpar(indalfap)*rpar(inddeltp)
c
c... compute current viscous strain rate
c
@@ -684,7 +669,7 @@
c
call dcopy(nstr,rpar(indconst),ione,r,ione)
call daxpy(nstr,rmtdt,betatdt,ione,r,ione)
- call dspmv("u",nstr,sub,rpar(inddmati),scur,ione,one,r,ione)
+ call dspmv("u",nstr,one,rpar(inddmati),scur,ione,one,r,ione)
c
return
end
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-15 17:15:26 UTC (rev 6261)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/newt.f 2007-03-15 19:38:25 UTC (rev 6262)
@@ -4,6 +4,9 @@
c... routine to find zero of a function using Newton's method with
c line search.
c Adapted from Numerical Recipes.
+c** Note that in this version, the Jacobian is assumed to always
+c be symmetric and stored in packed triangular form, with dimensions
+c of 6x6.
c
include "implicit.inc"
c
@@ -11,6 +14,7 @@
c
include "nconsts.inc"
include "rconsts.inc"
+ include "ndimens.inc"
integer NP,MAXITS
double precision TOLF,TOLMIN,TOLX,STPMX
parameter(NP=6,MAXITS=100,TOLF=1.0d-8,TOLMIN=1.0d-10,TOLX=1.0d-12,
@@ -34,7 +38,7 @@
double precision fvec(NP)
CU USES fdjac,fmin,lnsrch,lubksb,ludcmp
INTEGER i,its,j,indx(NP)
- double precision d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP)
+ double precision d,den,f,fold,stpmax,sum,temp,test,fjac(nddmat)
double precision g(NP),p(NP),xold(NP),fmin
external fmin
c
@@ -42,38 +46,27 @@
if(ierr.ne.izero) return
test=zero
do i=1,n
- if(abs(fvec(i)).gt.test)test=abs(fvec(i))
+ test=max(test,abs(fvec(i)))
end do
if(test.lt.0.01d0*TOLF)return
sum=dnrm2(n,x,ione)
stpmax=STPMX*max(sum,dble(n))
do its=1,MAXITS
- call jaccmp(x,n,fvec,NP,fjac,rpar,nrpar,ipar,nipar,
+ call jaccmp(x,n,fvec,nddmat,fjac,rpar,nrpar,ipar,nipar,
& ierr,errstrng)
-c**** I should replace this with a BLAS call
- do i=1,n
- sum=zero
- do j=1,n
- sum=sum+fjac(j,i)*fvec(j)
- end do
- g(i)=sum
- end do
+ call dspmv("u",n,one,fjac,fvec,ione,zero,g,ione)
call dcopy(n,x,ione,xold,ione)
fold=f
do i=1,n
p(i)=-fvec(i)
end do
-c**** these can be replaced by LAPACK calls -- eventually, it will
-c**** probably make more sense to use symmetric packed format for
-c**** all matrices.
- call ludcmp(fjac,n,NP,indx,d)
- call lubksb(fjac,n,NP,indx,p)
- call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin,rpar,nrpar,
- & ipar,nipar,funcv,ierr,errstrng)
+ call dppsv('u',n,1,fjac,p,n,ierr)
+ call lnsrch(n,xold,fold,g,p,x,f,fvec,stpmax,check,fmin,
+ & rpar,nrpar,ipar,nipar,funcv,ierr,errstrng)
if(ierr.ne.izero) return
test=zero
do i=1,n
- if(abs(fvec(i)).gt.test)test=abs(fvec(i))
+ test=max(test,abs(fvec(i)))
end do
if(test.lt.TOLF)then
check=.false.
@@ -84,7 +77,7 @@
den=max(f,half*dble(n))
do i=1,n
temp=abs(g(i))*max(abs(x(i)),one)/den
- if(temp.gt.test)test=temp
+ test=max(test,temp)
end do
if(test.lt.TOLMIN)then
check=.true.
@@ -96,10 +89,11 @@
test=zero
do i=1,n
temp=(abs(x(i)-xold(i)))/max(abs(x(i)),one)
- if(temp.gt.test)test=temp
+ test=max(test,temp)
end do
if(test.lt.TOLX)return
end do
-c**** replace this with an error code
- pause 'MAXITS exceeded in newt'
+ ierr=120
+ errstrng="newt"
+ return
END
More information about the cig-commits
mailing list