c PROGRAM VC_SG_COMPRESS.F c c Computation of Stress Green's functions (influence functions) c for interaction of fault segments. Present code c allows possibility of both elastic and viscoelastic c Greens functions according to the metod of c Rundle (1988a). This code is run on the output c file from VC_STRESS_GREEN.F and results in a c condensation of the m-dependent Laplace-transform c image terms in the viscoelastic stress Greens c function. c c Modified July 5, 2002 c c ************************************************************* c c General Description: c c This is the code that computes the stress green's functions c according to the methods described in Rundle (1988a,b) for c a prescribed group of vertical strike slip fault segments c embedded in an elastic half space, or in an elastic layer c overlying an viscoelastic half space. The stress greens c function specifies the stress change on fault j due to c a unit of slip on fault i. Relative geometry of the fault c segments, as well as elastic or viscoelastic properties of c the medium, are important. The output file of this code is c used as input file to VC_SG_Compress.f, which compresses c the viscoelastic part into a sum of 8 decaying exponentials c (collocation technique) for each fault segment's influence c functions. The output file of VC_Stress_Adjust.f is then c used as input to the first Virtual California code, c VIRSCA_INIT_SER.f. The output of that code is then used as c input to VIRSCA_SER.f c c c References: c c The basic ideas behind these codes are described in the c following series of papers: c c Rundle, JB, A physical model for earthquakes, 1., Fluctuations c and interactions, J. Geophys. Res., 93, 6237 (1988a) c c Rundle, JB, A physical model for earthquakes, 2. Application c to Southern California, J. Geophys. Res., 93, 6255 c (1988b) c c Rundle, JB, Linear pattern dynamics in nonlinear threshold c systems, Phys. Rev. E, 61, 2418 (2001) c c Rundle, PB, JB Rundle, KF Tiampo, JSS Martins, S McGinnis, and c W Klein, Phys. Rev. Lett., 87, 148501 (2001) c c Rundle, JB, PB Rundle, W Klein, JSS Martins, KF Tiampo, c A Donnellan and LH Kellogg, GEM plate boundary c simulations for the plate boundary observatory: c Understanding the physics of earthquakes c on complex fault systems, Pure. Appl. Geophys., c in press (2002) c c ************************************************************* c c Virtual California Code Physics: c c Virtual California uses topologically realistic networks of c independent fault segments that are mediated by elastic c interactions. (Note that earlier versions had viscoelastic c interactions as well, but such are not included in the present c code). VC is a "backslip" model, inasmuch as the plate c tectonic stress increases is produced by means of applying a c negative ("backslip") velocity to each segment whose magnitude c is that of the long-term rate of slip on the segment. Since c "positive slip" reduces the stress on a fault segment, c "negative slip" due to the backslip increases the stress. c On each time step, all faults are checked to determine whether c the shear stress has reached the failure threshold. Once c at least one segment reaches the threshold, the "long time c steps" stop, and "short (failure) time steps" (a.k.a. Monte c Carlo Sweeps, or mcs) begins. An mcs begins with a check of c of each site to determine whether it has failed, followed by c a parallel updating of each segment. An update of a segment c consists of increasing the sudden seismic slip on each segment c so that the stress of the segment, considered in isolation, c drops to a residual value, plus or minus a random overshoot/ c undershoot. The elastic stress on all segments is then c recalculated, and another mcs is carried out. This iterative c process repeats until all segments are below the failure c threshold, at which time the mcs time steps cease and the c long plate tectonic time steps begin again. Note that VC c also includes a stress-dependent "precursory slip", or c "stress leakage" of the type that has been observed in c laboratory experiments by TE Tullis (1996) and S Karner c and C Marone (2001). The physics of this process is that c as the stress on a segment increases, a small amount of stable c sliding occurs that is proportional to the level of the stress c above the residual. Lab experiments and field data suggest c that the parameter called "alpha" (in Rundle et al. 2001) c is of the order of a few percent. This parameter is called c facsdp() or facrat() or facsd() in this code. c c In addition, as described in Rundle (1988), the fault system c topology + elasticity + failure law may lead to an unstable, c "runaway" dynamics due to the presence of positive c eigenvalues, so we introduce a nonlinear, "self-adapting" or c "self-correcting" term. The VC initializer code estimates c the magnitude of these terms for each segment using a simple c but unrealistic dynamics. They are then readjusted by c running the VC code with real dynamics until all startup c transients (unrealistically large stresses) disappear. The c time scale for the transients to disappear depends on the c complexity of the problem. Physically, one can imagine c that the system evolves on a rough energy landscape, and c on average sits at the bottom of a free energy well. Large c earthquakes may tend to displace the system from the well c allowing non-stationary evolution to occur before the system c settles into a new well. The configuration of the well is c determined by the magnitude of the positive eigenvalues c discussed above. The eigenvalues for each segment are c proportional to the derivative of the shear stress with c respect to the slip on each fault segment. c c In this code, VC_SG_Compress.f, the m-dependent terms c which depend on Laplace transform variable s, are fit c to a series of 8 decaying exponentials in order to condense c the viscoelastic stress Green's functions further. We use c the method discussed, for example, in RA Schapery (Proc. c 4th US Nat. Cong. Appl. Mech., 1075, 1962), and discussed c in appendix C of the book by: RM Christensen, Theory of c Viscoelasticity, An Introduction, Academic Press, New York, c 1971. Basically, the method utilizes a matrix inversion c to find best-fitting coefficients to approximate the m- c dependent terms. c c ************************************************************* c c Program Modules: c c This code uses one main program and one subroutine c module, as well as a number of subroutine modules c from the LINPACK library written by Jack Dongarra, c and available from http://www.netlib.org/linpack/index.html c These LINPACK routines perform the matrix inversions. c c These modules are: c c Main code -- Reads in output file from c VC_Stress_Green.f, and calls c subroutine TREVNI() which inverts c the m-dependent series to find c the coefficients. This part of c code also defines the set of linear c equations to be inverted. c c Subroutine TREVNI() -- Subroutine that calls the LINPACK c routines to solve Ax=b. c c Various LINPACK routines that are used to solve the set of c linear equations Ax=b. c c ************************************************************* c c Dimensions from INCLUDE statement: c c NF -- MAX number of faults -- should be same as in EQPARM.FOR c file c NTSTP -- Max number of time steps c NFT -- Same as NTSTP c NF1 -- Same as NF c NF2 -- Same as NF c NFF -- NF * NF c NMF -- > NF * NF * 8 c NPT -- > NF * 36 c c c ************************************************************* c c Input & Output Files: c c c 1. INPUT FILE. The primary input (other than program c interactive inputs from screen) is a file that c contains the shear stress and normal stress Green's c functions, or influence coefficients that have c been computed from VC_STRESS_GREEN.F Both elastic c and viscoelastic Green's functions are stored. c Convention is "+" is right lateral shear stress, c "-" is left lateral shear stress. c c The primary input variables are: c c fname5 (fndum) - Character string name of input file c numflt (nfault) - Number of fault segments c mend - Number of viscoelastic images c hplate - Elastic plate thickness in km c vplatx - x-coordinate of far field plate tectonic c velocity in cm/yr c vplaty - y-coordinate of far field plate tectonic c velocity in cm/yr c bulkm - Elastic bulk modulus in bars c c rhofcc(NF) - Pressure at center of fault segment in c bars c c Ssh(NFF)- Shear stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Snorm(NFF)- Normal stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Sshm(NMF)- Shear stress viscoelastic Green's function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c Snormm(NMF)- Normal stress viscoelastic Green function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c area(NF) - areas of segments in km^2 c c svec(NF) - represents the total shear stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c snvec(NF) - Represents the total normal stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c slpvel(NF) - Slip velocity in cm/yr. Convention is c "+" (>0) for right lateral strike slip, "-" c (<0) for left lateral strike slip. c c vcovec(NFF) - Used in computing svec(NF). Represents c the total influence of a unit of slip on c fault segment j on the shear stress change c on fault segment i. c c vnvec(NFF) - Used in computing snvec(NF). Represents c the total influence of a unit of slip on c fault segment j on the normal stress change c on fault segment i. c c c 2. OUTPUT FILE. The primary output (other than program c diagnostic printouts to screen) is a file that c contains the shear stress and normal stress Green's c functions, or influence coefficients that have c been computed from VC_Stress_Green.f Both elastic c and viscoelastic Green's functions are stored, but c the viscoelastic Green's functions have now been c compressed to 8 terms for each segment. c Convention is "+" is right lateral shear stress, c "-" is left lateral shear stress. c c The primary outputs variables are: c c fname5 - Character string name of input file c numflt - Number of fault segments c mend - Number of viscoelastic terms = 8 now c hplate - Elastic plate thickness in km c vplatx - x-coordinate of far field plate tectonic c velocity in cm/yr c vplaty - y-coordinate of far field plate tectonic c velocity in cm/yr c bulkm - Elastic bulk modulus in bars c c c rhofcc(NF) - Pressure at center of fault segment in c bars c timpol(8) - Time constants for the 8 terms in the c viscoelastic sum of decaying exponentials for c each fault segment. Units of yrs. c Ssh(NFF)- Shear stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Snorm(NFF)- Normal stress elastic Green's function for c influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j. c c Sshm(NMF)- Shear stress viscoelastic Green's function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c Snormm(NMF)- Normal stress viscoelastic Green function c for influence of segment i on segment j. Units c are in inverse cm (cm^-1). Will be multiplied c by amu prior to using in Virtual California. c Note that the N*N*mend coefficients are stored c by columns as usual, with first column being c source fault i, all N target faults j, all c mend images. c c area(NF) - areas of segments in km^2 c c svec(NF) - represents the total shear stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c snvec(NF) - Represents the total normal stress change c per year, in inverse yr, (yr^-1) c (i.e., normalized by amu in bars), due to c the loading velocity on all the faults. c c slpvel(NF) - Slip velocity in cm/yr. Convention is c "+" (>0) for right lateral strike slip, "-" c (<0) for left lateral strike slip. c c vcovec(NFF) - Used in computing svec(NF). Represents c the total influence of a unit of slip on c fault segment j on the shear stress change c on fault segment i. c c vnvec(NFF) - Used in computing snvec(NF). Represents c the total influence of a unit of slip on c fault segment j on the normal stress change c on fault segment i. c c ************************************************************* INCLUDE 'EQPARM_4.FOR' implicit real*8(a-h,o-z) character*20 fndum,fnstro, fnstrn dimension Sxy(NFF),Sn(NFF),Sxym(NMF), . Snm(NMF),timpol(8),Sxyw(NMF),Snw(NMF) dimension amat(8,8), svecxy(8), svecn(8) dimension ansvxy(8),ansvn(8), bmat(8,8) dimension gamma(NF), rhofac(NF), Sxyvec(NF), . Snvec(NF), Slpvel(NF), vcovec(NFF), vnvec(NFF) c c Read the fault coefficients c print *,' ' print *,' Enter the file name for the output' print *,' from program STRESS_COEFFICIENTS' read(*,'(a20)') fnstro c print *,' ' print *,' Enter the file name for OUTPUT' print *,' from this program' read(*,'(a20)') fnstrn c pole1=.25d0 pmult=2.d0 timpol(1)=pole1 do i=2,8 timpol(i)=timpol(i-1)*pmult end do c open(unit=10,name=fnstro,status='old') read(10,'(a20)') fndum read(10,*) nfault,mend,hplate,vplatx,vplaty,bulkm read(10,*) (rhofac(n), n=1,nfault) c c Read the fault coefficients c imax=nfault*nfault imaxm=nfault*nfault*mend read(10,*) (Sxy(n), n=1,imax) read(10,*) (Sn(n), n=1,imax) read(10,*) (Sxym(n), n=1,imaxm) read(10,*) (Snm(n), n=1,imaxm) read(10,*) (gamma(n), n=1,nfault) read(10,*) (Sxyvec(n), n=1,nfault) read(10,*) (Snvec(n), n=1,nfault) read(10,*) (Slpvel(n), n=1,nfault) read(10,*) (vcovec(n), n=1,imax) read(10,*) (vnvec(n), n=1,imax) close(10) c c Take only first ten terms because of gravity effect c mindx=mend itaska=1 itaskb=1 index=0 do i=1,nfault Sxyvec(i)=0.d0 Snvec(i)=0.d0 end do do ic=1,nfault print *,' CALCULATIONS FOR PATCH --',ic do ir=1,nfault ind=(ic-1)*nfault + ir indm=(ic-1)*mend*nfault + (ir-1)*mend ind8=(ic-1)*8*nfault + (ir-1)*8 c do i=1,8 si=1.d0/timpol(i) do j=1,8 sj=1.d0/timpol(j) amat(i,j)=sj/(si*(si+sj)) bmat(i,j)=amat(i,j) end do fact=1.d0 Svecxy(i)=0.d0 Svecn(i)=0.d0 do m=1,mindx indmx=indm + m fact=(1.d0/(1.d0+si))*fact if(dabs(fact) .le. 1.d-20) fact=0.d0 svecxy(i)=fact*Sxym(indmx) + svecxy(i) svecn(i)=fact*Snm(indmx) + svecn(i) end do Svecxy(i)=Svecxy(i)/si Svecn(i)=Svecn(i)/si end do c call TREVNI(amat,svecxy,ansvxy,itaska) call TREVNI(bmat,svecn,ansvn,itaskb) c do m=1,8 ind8x=ind8 + m Sxyw(ind8x)=ansvxy(m) Snw(ind8x)=0.d0 Snw(ind8x)=ansvn(m) end do c axy=Sxy(ind) an=Sn(ind) do i=1,8 axy=axy + ansvxy(i) an=an + ansvn(i) end do index=index + 1 vcovec(index)=axy vnvec(index)=an Sxyvec(ir)=Sxyvec(ir) - axy*slpvel(ic) Snvec(ir)=Snvec(ir) - an*slpvel(ic) end do end do c c Write the file c mend=8 open(unit=20,name=fnstrn,status='new') write(20,'(a20)') fndum write(20,*) nfault,mend,hplate,vplatx,vplaty,bulkm write(20,*) (timpol(n), n=1,8) write(20,*) (rhofac(n), n=1,nfault) c c Write the fault coefficients c imax=nfault*nfault imaxm=nfault*nfault*8 write(20,*) (Sxy(n), n=1,imax) write(20,*) (Sn(n), n=1,imax) write(20,*) (Sxyw(n), n=1,imaxm) write(20,*) (Snw(n), n=1,imaxm) write(20,*) (gamma(n), n=1,nfault) write(20,*) (Sxyvec(n), n=1,nfault) write(20,*) (Snvec(n), n=1,nfault) write(20,*) (Slpvel(n), n=1,nfault) write(20,*) (vcovec(n), n=1,imax) write(20,*) (vnvec(n), n=1,imax) close(20) c stop end c subroutine TREVNI(Amat,sdrop,delta,itask) c c ************************************************************* c c c INTERFACE: c Amat(i,j) - The matrix A for each fault segment c sdrop(i) - The right hand side, representing c parameter b, for each fault segment. c delta(i) - The best fitting coefficient to be c determined, representing parameter x, for c each fault segment c itask - An error code generated by LINPACK c c RETURNED VALUE: c c DESCRIPTION: c This subroutine computes the best-fitting coefficients c for the m-dependent image terms for each fault segment, c using math. routines c from the LINPACK library, available c http://www.netlib.org/linpack/index.html See c the subroutine source code there, or below for details c of numerical methods. c c Method: For each fault segment, we take the mend terms c from Sxym and Snm, and fit these to 8 terms by solving an c equation Ax = b. c c ************************************************************* c implicit real*8(a-h,o-z) dimension Amat(8,8),sdrop(8),delta(8) dimension work(8),iwork(8) c c Invert the system of equations c ijob = 0 call dgeco(Amat,8,8,iwork,rcond,work) call dgesl(Amat,8,8,iwork,sdrop,ijob) c if(ind .lt. 0) then print *,' ** IND in DGEFS .lt. 1 **' print *,' IND =', IND end if c do i=1,8 delta(i)=sdrop(i) end do c return end c c c ************************************************************* double precision function dasum(n,dx,incx) c ------------------------------------------------------------------c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c ------------------------------------------------------------------c double precision dx(1),dtemp integer i,incx,m,mp1,n,nincx dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c --------------------------------- c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c clean-up loop c ----------------------------- c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end c c ************************************************************* c subroutine daxpy(n,da,dx,incx,dy,incy) c---------------------------------------------------------------c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c---------------------------------------------------------------c double precision dx(1),dy(1),da integer i,incx,incy,m,mp1,n if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ----------------------------------------------- c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c clean-up loop c ----------------------------------- c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end c c ************************************************************* c double precision function ddot(n,dx,incx,dy,incy) c-------------------------------------------------------------c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c-------------------------------------------------------------c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ----------------------------------------------- c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c clean-up loop c ----------------------------------- c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end c c ************************************************************* c subroutine dgeco(a,lda,n,ipvt,rcond,z) integer lda,n,ipvt(1) double precision a(lda,1),z(1) double precision rcond c c INTERFACE: c a double precision(lda, n) c the matrix to be factored. c lda integer c the leading dimension of the array a . c n integer c the order of the matrix a . c c RETURN VALUE: c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c ipvt integer(n) c an integer vector of pivot indices. c rcond double precision c an estimate of the reciprocal condition of a . c for the system a*x = b , relative perturbations c in a and b of size epsilon may cause c relative perturbations in x of size epsilon/rcond . c if rcond is so small that the logical expression c 1.0 + rcond .eq. 1.0 c is true, then a may be singular to working c precision. in particular, rcond is zero if c exact singularity is detected or the estimate c underflows. c z double precision(n) c a work vector whose contents are usually unimportant. c if a is close to a singular matrix, then z is c an approximate null vector in the sense that c norm(a*z) = rcond*norm(a)*norm(z) . c c DESCRIPTION: c dgeco factors a double precision matrix by gaussian elimination c and estimates the condition of the matrix. c c if rcond is not needed, dgefa is slightly faster. c to solve a*x = b , follow dgeco by dgesl. c to compute inverse(a)*c , follow dgeco by dgesl. c to compute determinant(a) , follow dgeco by dgedi. c to compute inverse(a) , follow dgeco by dgedi. c c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c linpack dgefa c blas daxpy,ddot,dscal,dasum c fortran dabs,dmax1,dsign c c internal variables c double precision ddot,ek,t,wk,wkm double precision anorm,s,dasum,sm,ynorm integer info,j,k,kb,kp1,l c c c compute 1-norm of a c anorm = 0.0d0 do 10 j = 1, n anorm = dmax1(anorm,dasum(n,a(1,j),1)) 10 continue c c factor c call dgefa(a,lda,n,ipvt,info) c c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . c trans(a) is the transpose of a . the components of e are c chosen to cause maximum local growth in the elements of w where c trans(u)*w = e . the vectors are frequently rescaled to avoid c overflow. c c solve trans(u)*w = e c ek = 1.0d0 do 20 j = 1, n z(j) = 0.0d0 20 continue do 100 k = 1, n if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k)) if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30 s = dabs(a(k,k))/dabs(ek-z(k)) call dscal(n,s,z,1) ek = s*ek 30 continue wk = ek - z(k) wkm = -ek - z(k) s = dabs(wk) sm = dabs(wkm) if (a(k,k) .eq. 0.0d0) go to 40 wk = wk/a(k,k) wkm = wkm/a(k,k) go to 50 40 continue wk = 1.0d0 wkm = 1.0d0 50 continue kp1 = k + 1 if (kp1 .gt. n) go to 90 do 60 j = kp1, n sm = sm + dabs(z(j)+wkm*a(k,j)) z(j) = z(j) + wk*a(k,j) s = s + dabs(z(j)) 60 continue if (s .ge. sm) go to 80 t = wkm - wk wk = wkm do 70 j = kp1, n z(j) = z(j) + t*a(k,j) 70 continue 80 continue 90 continue z(k) = wk 100 continue s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) c c solve trans(l)*y = w c do 120 kb = 1, n k = n + 1 - kb if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1) if (dabs(z(k)) .le. 1.0d0) go to 110 s = 1.0d0/dabs(z(k)) call dscal(n,s,z,1) 110 continue l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t 120 continue s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) c ynorm = 1.0d0 c c solve l*v = y c do 140 k = 1, n l = ipvt(k) t = z(l) z(l) = z(k) z(k) = t if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1) if (dabs(z(k)) .le. 1.0d0) go to 130 s = 1.0d0/dabs(z(k)) call dscal(n,s,z,1) ynorm = s*ynorm 130 continue 140 continue s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) ynorm = s*ynorm c c solve u*z = v c do 160 kb = 1, n k = n + 1 - kb if (dabs(z(k)) .le. dabs(a(k,k))) go to 150 s = dabs(a(k,k))/dabs(z(k)) call dscal(n,s,z,1) ynorm = s*ynorm 150 continue if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k) if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0 t = -z(k) call daxpy(k-1,t,a(1,k),1,z(1),1) 160 continue c make znorm = 1.0 s = 1.0d0/dasum(n,z,1) call dscal(n,s,z,1) ynorm = s*ynorm c if (anorm .ne. 0.0d0) rcond = ynorm/anorm if (anorm .eq. 0.0d0) rcond = 0.0d0 return end c c ************************************************************* c subroutine dgedi(a,lda,n,ipvt,det,work,job) integer lda,n,ipvt(1),job double precision a(lda,1),det(2),work(1) c c INTERFACE: c a double precision(lda, n) c the output from dgeco or dgefa. c lda integer c the leading dimension of the array a . c n integer c the order of the matrix a . c ipvt integer(n) c the pivot vector from dgeco or dgefa. c work double precision(n) c work vector. contents destroyed. c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c RETURN VALUE: c a inverse of original matrix if requested. c otherwise unchanged. c det double precision(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. dabs(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c ERROR CONDITION: c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if dgeco has set rcond .gt. 0.0 or dgefa has set c info .eq. 0 . c c DESCRIPTION: c dgedi computes the determinant and inverse of a matrix c using the factors computed by dgeco or dgefa. C c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,dscal,dswap c fortran dabs,mod c c internal variables c double precision t double precision ten integer i,j,k,kb,kp1,l,nm1 c c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0d0 det(2) = 0.0d0 ten = 10.0d0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (det(1) .eq. 0.0d0) go to 60 10 if (dabs(det(1)) .ge. 1.0d0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0d0 go to 10 20 continue 30 if (dabs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0d0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n a(k,k) = 1.0d0/a(k,k) t = -a(k,k) call dscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0.0d0 call daxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = 0.0d0 110 continue do 120 j = kp1, n t = work(j) call daxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end c c ************************************************************* c subroutine dgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info double precision a(lda,1) c c c INTERFACE: c a double precision(lda, n) c the matrix to be factored. c lda integer c the leading dimension of the array a . c n integer c the order of the matrix a . c c RETURN VALUE: c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c ipvt integer(n) c an integer vector of pivot indices. c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that dgesl or dgedi will divide by zero c if called. use rcond in dgeco for a reliable c indication of singularity. c c DESCRIPTION: c dgefa factors a double precision matrix by gaussian elimination. c c dgefa is usually called by dgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for dgeco) = (1 + 9/n)*(time for dgefa) . c c REFERENCES: c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,dscal,idamax c c internal variables c double precision t integer idamax,j,k,kp1,l,nm1 c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = idamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0d0/a(k,k) call dscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0d0) info = n return end c c ************************************************************* c subroutine dgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job double precision a(lda,1),b(1) c c c INTERFACE: c a double precision(lda, n) c the output from dgeco or dgefa. c lda integer c the leading dimension of the array a . c n integer c the order of the matrix a . c ipvt integer(n) c the pivot vector from dgeco or dgefa. c b double precision(n) c the right hand side vector. c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c RETURN VALUE: c b the solution vector x . c c ERROR CONDITION: c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgeco has set rcond .gt. 0.0 c or dgefa has set info .eq. 0 . c c c DESCRIPTION: c dgesl solves the double precision system c a * x = b or trans(a) * x = b c using the factors computed by dgeco or dgefa. c c to compute inverse(a) * c where c is a matrix c with p columns c call dgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c REFERENCES: c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,ddot c c internal variables c double precision ddot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end c c ************************************************************* c subroutine dscal(n,da,dx,incx) c----------------------------------------------------------------c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c----------------------------------------------------------------c double precision da,dx(1) integer i,incx,m,mp1,n,nincx if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c --------------------------------- c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c clean-up loop c ----------------------------- c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end c c ************************************************************* c subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end c c ************************************************************* c integer function idamax(n,dx,incx) c-----------------------------------------------------------------c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c-----------------------------------------------------------------c double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c --------------------------------- c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c ----------------------------- c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end