c c ***************************************************************** c * * c * PROC * c * program for processing sensitivity matrices * c * by principal component analysis * c * * c ***************************************************************** c c PROC is a member of KINAL , a program package for the c analysis of kinetic reaction mechanisms c c The following subroutines belong to PROC main program: c c EQS SDIAG2 ORDER FUN JFY JFP QSA c SOL DEC RATE RTIME c ----------------------------------------------------------------- c c mode 0 S matrix investigation c 1 F matrix investigation c 2 QS matrix investigation c c ny - number of species ( max 50 ) | in the c nr - number of reactions ( max 90 ) | original c np - number of parameters ( max 90 ) | data c nt - number of time points ( max 60 ) | file c c nyo - number of "observed" concentrations c npo - the first npo parameters in the data file are "observed" c the parameters are assumed to be in increasing order c nto - number of "observed" time points c lis - if lis.eq.0 listing of reactions is forbidden c c iny - serial numbers of "observed" species c int - serial numbers of "observed" time points in increasing order c i refers to the ith data block or to the (i+1)th time point c c formula - identification of species c y - concentrations c p - parameters c ta - time points ( ta(1) is the initial time ) c plam - lambda parameter of QSA computing c c ndp - leading dimension of matrix s c c ----------------------------------------------------------------- c c npp = np * (np+1) / 2 c c character*8 formula(ny) c dimension ta(nt),y(ny),dy(ny),yx(ny),int(nt),iny(ny),s(np,np),sTs(npp) c dimension ip(np),ipp(np),b(np),bb(np),ifn(np),ifon(np,25),fy(ny,ny) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c c ----------------------------------------------------------------- c c written by Tamas Turanyi c Department of Physical Chemistry c E”tv”s University c c Address: H-1518 Budapest-112, P.O.Box 32, Hungary c Phone : (36-1) 209-0555 ext. 1109 c Fax : (36-1) 209-0602 c E-mail : turanyi@garfield.chem.elte.hu c WWW : www-PhCh.chem.elte.hu c c ----------------------------------------------------------------- c c refer to: T. Turanyi, Comp.Chem., Vol. 14, pp. 253-254, 1990 c c ----------------------------------------------------------------- c implicit real*8(a-h,o-z) implicit integer*2 (i-n) character*1 c character*3 list character*8 formula(50),fal,comp character*10 datum,times character*14 fin,fout character*25 chm(3) character*80 rem dimension ta(60),y(50),dy(50),yx(50),int(60),iny(50) dimension ip(90),ipp(90),b(90),bb(90),ifn(90) common /w1/ fy(50,50),sTs(4095),ifon(90,25) common /w2/ s(90,90) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c ndp= 90 nca= 200 ncb= 300 npp = ndp * (ndp+1) / 2 c comp= 'data No.' chm(1)= ' S matrix investigation' chm(2)= ' F matrix investigation' chm(3)= ' QS matrix investigation' c write(*,*) 'Name of data file ?' read(*,'(a)') fin write(*,*) 'Name of output file ?' read(*,'(a)') fout c open(5,file=fin,status='old') open(6,file=fout,status='new') c read(5,103) rem read(5,100) ny,nr,np,nt,nyo,npo,nto,mode,lis read(5,100) (int(i),i=1,nto) read(5,100) (iny(i),i=1,nyo) read(5,101) (formula(i),y(i),i=1,ny) read(5,102) (p(i),i=1,nr) read(5,102) plam read(5,102) (ta(i),i=1,nt) c list=' no' if (lis.ne.0) list='yes' c if (mode.lt.0.or.mode.gt.2) then write(6,229) mode stop endif c write(6,200) rem,chm(mode+1) write(6,201) ny,nr,np,nt,nyo,npo,nto,list if(mode.eq.2) write(6,231) plam c if (nyo.gt.ny) then write(6,218) stop endif if (npo.gt.np) then write(6,219) stop endif if (nto.gt.(nt-1)) then write(6,220) stop endif c do 1 i=1,nto if(int(i).lt.1.or.int(i).gt.(nt-1)) then write(6,202) int(i) stop endif if (i.eq.1) goto 1 if (int(i).le.int(i-1)) then write(6,222) stop endif 1 continue do 2 i=1,nyo if(iny(i).lt.1.or.iny(i).gt.ny) then write(6,203) iny(i) stop endif 2 continue c write(6,204) (i,p(i),i=1,nr) write(6,205) (i,ta(int(i)+1),i=1,nto) write(6,206) (formula(iny(i)),i=1,nyo) c nisa= 0 nisb= 0 3 read(5,104) k,l,ns if(k.eq.0) goto 4 if(k.gt.nr.or.l.gt.ny.or.ns.gt.3) write(6,207) k,l,ns if(k.le.0.or.l.le.0) write(6,207) k,l,ns if(ns.eq.0) ns = 1 nisa= nisa + 1 if (nisa.gt.nca) then write(6,227) nca stop endif isa(1,nisa) = k isa(2,nisa) = l isa(3,nisa) = ns goto 3 4 read(5,105) k,l,sn if(k.eq.0) goto 5 if(k.gt.nr.or.l.gt.ny.or.sn.gt.3.) write(6,208) k,l,sn if(k.le.0.or.l.le.0) write(6,208) k,l,sn nisb = nisb + 1 if (nisb.gt.ncb) then write(6,228) ncb stop endif isb(1,nisb) = k isb(2,nisb) = l sb(nisb) = sn goto 4 5 continue c if (lis.ne.0.and.nr.gt.7) write(6,209) if (lis.ne.0) call eqs(formula) c call rtime(datum,times) write(6,210) datum,times c do 19 i=1,npp 19 sTs(i)= 0. do 31 i=1,ndp 31 ip(i)= i c do 11 it=1,nto 6 read (5,101) fal if (fal.eq.comp) goto 7 goto 6 7 read (5,106) ic if(ic.eq.int(it)) goto 8 goto 6 c 8 continue read (5,107) t,(y(i),i=1,ny) c call fun(y,dy) write(6,211) int(it),t write(6,217)(formula(i), y(i),i=1,ny) write(6,230)(formula(i),dy(i),i=1,ny) c do 35 i=1,ny 35 yx(i)= y(i) c if(mode-1) 27,28,29 27 do 9 j=1,npo read (5,107) (s(i,j),i=1,ny) 9 continue goto 30 28 call jfp(nr,ndp,y,s,ip) do 16 i=1,ny 16 yx(i)= dy(i) goto 30 29 call qsa(ny,nr,ndp,y,s,fy,plam,ip) c 30 continue c iv= 0 do 10 i=1,npo do 10 j=1,i iv= iv+1 do 10 k=1,nyo l=iny(k) if (dabs(yx(l)).lt.1.d-50) goto 10 sTs(iv)= sTs(iv) + s(l,i)*s(l,j)/yx(l)/yx(l)*p(i)*p(j) 10 continue 11 continue iv=0 do 15 i=1,npo do 15 j=1,i iv=iv+1 15 s(i,j)= sTs(iv) c c hereafter sTs is a working array c c overall sensitivities and reaction rate c c do 32 i=1,nr 32 sTs(i)= 0. do 33 i=1,npo 33 sTs(i)= s(i,i) call order(npo,ip,sTs,b) call rate(y,sTs) call order(nr,ipp,sTs,bb) write(6,232) npr= max0(nr,npo) do 34 i=1,npr 34 write(6,233) i,ip(i),b(i),ipp(i),bb(i) c c evaluation of eigenvectors and eigenvalues by c the SDIAG2 routine c call sdiag2(ndp,npo,s,bb) c note : rank of sTs matrix.le.nsr nsr= min0(npo,nto*nyo) c c the principal components c write(6,214) do 17 i=1,nsr if(bb(i).lt.1.d-10) goto 12 write(6,212) i,bb(i) call order (npo,ip,s(1,i),b) do 14 j=1,npo if(dabs(b(j)).lt..05 ) goto 14 write(6,213) j,ip(j),b(j) 14 continue 17 continue c c choosing the thresholds for mechanism reduction c 12 continue do 20 i=1,ndp 20 ifn(i)= 0 do 21 i=1,ndp do 21 j=1,25 21 ifon(i,j)= 0 c write(*,225) read(*,*) vse write(*,226) read(*,*) vsv c do 25 i=1,nsr if(bb(i).lt.vse) goto 26 do 24 j=1,npo if(dabs(s(j,i)).lt.vsv) goto 24 ifn(j) = ifn(j) + 1 ifon(j,ifn(j)) = i 24 continue 25 continue 26 write(6,221) vse,vsv do 18 i=1,npo write(6,215) i if(ifn(i).eq.0) goto 18 ifl=ifn(i) write(6,216) (ifon(i,j),j=1,ifl) 18 continue c 23 write(*,224) read(*,223) c if(c.eq.'y') goto 12 if(c.eq.'n') goto 22 goto 23 c 22 call rtime(datum,times) write(6,210) datum,times write(6,209) c 100 format(16i5) 101 format(a8,2x,f10.2) 102 format (8f10.2) 103 format(a80) 104 format(10i5) 105 format(2i5,f5.2) 106 format (i3) 107 format (6x,6e12.5) 108 format (i3,3x,6e12.5/10(6x,6e12.5/)) c 200 format(1x,a80//1x,a25) 201 format(//,3x,' Original data file :'//3x, 1 'number of species : ',i3/3x, 2 'number of reactions : ',i3/3x, 3 'number of parameters : ',i3/3x, 4 'number of time points : ',i3///3x, 5 ' "Observed" data :'//3x, 6 'number of species : ',i3/3x, 7 'number of parameters : ',i3/3x, 8 'number of time points : ',i3//3x, 9 'list of reactions : ',a3/) 202 format (/2x,'Error : unrecognisable time point reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 203 format (/2x,'Error : unrecognisable obs. conc. reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 204 format (/2x,' Rate coefficients :'//100(4(2x,i5,1h.,1pe12.3)/)) 205 format (/2x,' Time points :'//100(4(2x,i5,1h.,1pe12.3)/)) 206 format (/2x,' Observed concentrations :'// 10(7(2x,a8)/)) 207 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i 3,' species :',i3,' stoich.coeff. :',i3) 208 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i3,' species :',i3,' stoich.coeff. :',f5.2) 209 format(1h1) 210 format(//20x,' Date : ',a10,' Time : ',a10//) 211 format (/' Reading data no.',i3, 1 ' Time : ',1pe15.3) 212 format(//3x,'No',i4,' eigenvalue :',1pe15.5/9x,' eigenvector :'/) 213 format(5x,i3,'.',i5,f6.3) 214 format(1h1,5x,'Principal component analysis :'//) 215 format(5x,i3) 216 format(1h+,10x,20i3,5(/11x,20i3)) 217 format (/2x,' Concentrations :'// 100(3(2x,a8,2h =,1pd13.5)/)) 218 format (/2x,'Error : nyo greater than ny'// 1 2x,' --- PROGRAM TERMINATED ---'///) 219 format (/2x,' Error : npo greater than np'// 1 2x,' --- PROGRAM TERMINATED ---'///) 220 format (/2x,'Error : nto greater than nt'// 1 2x,' --- PROGRAM TERMINATED ---'///) 221 format(1h1,' Proposal for mechanism reduction :'/ / 1 /2x,' Threshold value for eigenvalues :',1pe15.5, 2 /2x,' Threshold value for eigenvectors :',0pf12.5/) 222 format (/2x,'Error : time point references are not '/ 1 /2x,' in increasing order '/// 2 /2x,' --- PROGRAM TERMINATED ---'///) 223 format (1a1) 224 format (//' Another cut ? y / n') 225 format (//' Threshold value for eigenvalues ?') 226 format (/ ' Threshold value for eigenvectors ?') 227 format(//2x,'Error : redimension isa of ST common in main '// 1 2x,' and in the following subroutines '// 2 2x,' EQS , FUN , JFY , JFP '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '/ //) 228 format(//2x,'Error : redimension sb,isb of ST common in main'// 1 2x,' and in the following subroutines '// 2 2x,' EQS , FUN , JFY , JFP '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '///) 229 format(//2x,'Error : ',i5,' is not existing mode '// 1 2x,' --- PROGRAM TERMINATED --- '///) 230 format (/2x,' Production rates of species :'/ 1 /100(3(2x,a8,2h =,1pd13.5)/)) 231 format (/2x,' Lambda :',1pd13.5) 232 format (///5x,'No overall sens. reaction rates'//) 233 format (4x,i3,i6,1pe10.2,i8,1pe10.2) c close (6, status = 'keep' ) stop end c c *** KINAL ****************************************************** c * * c * EQS makes the list of reactions * c * * c **************************************************************** c c formula - IN: identification of species c subroutine eqs(formula) implicit integer*2 (i-n) implicit real*8 (a-h,o-z) character left(3) character*3 plus(10) character*8 blank,formula(1),aleft(3),aright(10) dimension right(10) c common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) data plus / 10*' + ' / c c the elements of isa ; isb,sb vectors are assumed c to be arranged in increasing order of reactions c iper = 6 blank =' ' jsa = 0 jsb = 0 do 1 k=1,nr il = 0 ir = 0 do 2 i=1,3 left(i) = ' ' 2 aleft(i)=blank do 3 i=1,10 right(i) = 0. 3 aright(i)=blank c 4 continue jsa= jsa + 1 if (jsa.gt.nisa) goto 7 if (isa(1,jsa).ne.k) goto 5 il= il + 1 if(il.gt.3) then il= 3 write(iper,100) k goto 7 endif aleft(il) = formula(isa(2,jsa)) if(isa(3,jsa).eq.2) left (il) = '2' if(isa(3,jsa).eq.3) left (il) = '3' goto 4 5 continue jsa = jsa - 1 7 continue jsb= jsb + 1 if (jsb.gt.nisb) goto 8 if (isb(1,jsb).ne.k) goto 6 ir= ir + 1 if(ir.gt.10) then ir= 10 write(iper,110) k goto 8 endif aright(ir) = formula(isb(2,jsb)) right(ir) = sb(jsb) goto 7 6 continue jsb = jsb - 1 8 continue c il1 = il - 1 ir3 = ir - 3 ir2 = min0(ir,2) ir1 = ir2 - 1 write(iper,120) k,(left(i),aleft(i),i=1,3), 1 (right(i),aright(i),i=1,ir2) write(iper,130) (plus(i),i=1,il1) write(iper,140) if(ir1.gt.0) write(iper,150) plus(1) if(ir .gt.2) then write(iper,160) (right(i),aright(i),i=3,ir) endif c 1 continue c 100 format(1x,'More then 3 species on the left side of eqs. No ',i3) 110 format(1x,'More then 10 species on the right side of eqs. No ',i3) 120 format(/1x,i3,'. ',3(a1,1x,a8,3x),2x,2(f4.2,1x,a8,3x)) 130 format(1h+,15x,3(a3,11x)) 140 format(1h+,41x,' --> ') 150 format(1h+,59x,a3,13x,a3) 160 format(2(/9x,4(' + ',f5.2,1x,a8))) 170 format(1h+,22x,8(a3,13x)) return end c c *** KINAL ****************************************************** c * * c * SDIAG2 computes eigenvectors and eigenvalues of a * c * symmetric matrix * c * * c **************************************************************** c c Method of QR transformations. c If the Euclidean norm of the rows varies STRONGLY most accurate c results may be obtained by permuting rows and columns to give c an arrangement with increasing norms of rows. c c Origin : Leibniz-Rechnenzentrum , Munich 1965 c c ---------------------------------------------------------------- c c x - IN: matrix to be diagonalized , its lower triangle has c to be given as ((x(i,j),j=1,i), i=1,n) c - OUT: matrix of eigenvectors c m - IN: leading dimension of matrix x c n - IN: dimension of matrix x c d - OUT: vector of eigenvalues c c ---------------------------------------------------------------- c c dimension e(n) c subroutine sdiag2(m,n,x,d) implicit real*8(a-h,o-z) implicit integer*2(i-n) dimension d(1),x(m,1),e(100) c c eps= minimum of all x such that 1+x is greater than 1 on the c computer eps= 1.0d-15 c c tol= inf / eps where inf = minimum of all positive x c representable within the computer c tol= 1.0d-300/1.0d-15 if (n.eq.1) goto 400 c c Householder's reduction c simulation of loop do 150 i=n,2,(-1) c do 150 ni=2,n ii= n+2-ni c c loop for recursive address calculation c do 150 i=ii,ii l= i-2 h= 0.0d0 g= x(i,i-1) if(l) 140,140,20 20 do 30 k=1,l 30 h= h+x(i,k)*x(i,k) s= h+g*g if (s.ge.tol) goto 50 h= 0.0d0 goto 140 50 if (h) 140,140,60 60 l= l+1 f= g g= dsqrt(s) if (f) 75,75,70 70 g= -g 75 h= s-f*g x(i,i-1)=f-g f= 0.0d0 do 110 j=1,l x(j,i)= x(i,j)/h s= 0.0d0 do 80 k=1,j 80 s= s+x(j,k)*x(i,k) j1= j+1 if (j1.gt.l) goto 100 do 90 k=j1,l 90 s= s+x(k,j)*x(i,k) 100 e(j)= s/h 110 f= f+s*x(j,i) f= f/(h+h) do 120 j=1,l 120 e(j)= e(j)-f*x(i,j) do 130 j=1,l f= x(i,j) s= e(j) do 130 k=1,j 130 x(j,k)= x(j,k)-f*e(k)-x(i,k)*s 140 d(i)= h 150 e(i-1)= g c c accumulation of transformation matrices c d(1)= x(1,1) x(1,1)= 1.0d0 do 220 i=2,n l= i-1 if (d(i)) 200,200,170 170 do 190 j=1,l s= 0.0d0 do 180 k=1,l 180 s= s+x(i,k)*x(k,j) do 190 k=1,l 190 x(k,j)= x(k,j)-s*x(k,i) 200 d(i)= x(i,i) x(i,i)= 1.0d0 do 220 j=1,l x(i,j)= 0.0d0 220 x(j,i)= 0.0d0 c c diagonalization of the tridiagonal matrix c b= 0.0d0 f= 0.0d0 e(n)= 0.0d0 do 340 l=1,n h= eps*(dabs(d(l))+dabs(e(l))) if (h.gt.b) b=h c c test for splitting c do 240 j=l,n if (dabs(e(j)).le.b) goto 250 240 continue c c test for convergence c 250 if (j.eq.l) goto 340 c c shift from upper 2*2 minor c 260 p= (d(l+1)-d(l))*0.5d0/e(l) r= dsqrt(p*p+1.0d0) p= p+dsign(r,p) h= d(l)-e(l)/p do 300 i=l,n 300 d(i)= d(i)-h f= f+h c c QR transformation c p= d(j) c= 1.0d0 s= 0.0d0 c c simulation of loop do 330 i=j-1,l,(-1) c j1= j-1 do 330 ni=l,j1 ii= l+j1-ni c c loop for recursive address calculation c do 330 i=ii,ii i1= i+1 g= c*e(i) h= c*p c c protection against underflow of exponents c if (dabs(p).lt.dabs(e(i))) goto 310 c= e(i)/p r= dsqrt(c*c+1.0d0) e(i1)= s*p*r s= c/r c= 1.0d0/r goto 320 310 c= p/e(i) r= dsqrt(c*c+1.0d0) e(i1)= s*e(i)*r s= 1.0d0/r c= c/r 320 p= c*d(i)-s*g d(i1)= h+s*(c*g+s*d(i)) do 330 k=1,n h= x(k,i1) x(k,i1)= x(k,i)*s+h*c 330 x(k,i)= x(k,i)*c-h*s e(l)= s*p d(l)= c*p if (dabs(e(l)).gt.b) goto 260 c c convergence c 340 d(l)= d(l)+f c c ordering of eigenvalues c ni= n-1 do 380 i=1,ni k= i p= d(i) j1= i+1 do 360 j= j1,n c if (d(j).le.p) goto 360 k= j p= d(j) 360 continue if (k.eq.i) goto 380 d(k)= d(i) d(i)= p do 370 j=1,n p= x(j,i) x(j,i)= x(j,k) 370 x(j,k)= p 380 continue return c c special treatment of case n= 1 c 400 d(1)= x(1,1) x(1,1)=1.0d0 return end c c *** KINAL ****************************************************** c * * c * ORDER puts into order the elements of a vector * c * according to their absolute value * c * * c **************************************************************** c c a - IN: vector of n elements (remains unchanged) c b - OUT: ordered vector c n - IN: dimension of vectors c ip - OUT: integer vector of n elements a(ip(i)) = b(i) c e - working area c c ---------------------------------------------------------------- c c dimension e(n) c subroutine order (n,ip,a,b) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension ip(1),a(1),b(1),e(100) do 1 i=1,n 1 e(i) = dsign(1.d+00,a(i)) c c rh must be greater then the greatest element of a rh = 1.d+50 do 2 i=1,n c c rm must be smaller then the smallest element of a rm = 1.d-50 do 3 j=1,n if(dabs(a(j)).ge.rm.and.dabs(a(j)).le.rh) goto 4 goto 3 4 if(i.eq.1) goto 5 do 6 k=2,i if(j.eq.ip(k-1)) goto 3 6 continue 5 rm=dabs(a(j)) im=j 3 continue b(i)=rm * e(im) rh=rm ip(i) = im 2 continue return end c c *** KINAL ****************************************************** c * * c * FUN computes the right hand side of kinetic diff. eq. * c * * c **************************************************************** c c y - IN: actual concentrations c dy - OUT: derivative of variables with respect to time c dy(i) = d y(i) / d t c r - working area c c ---------------------------------------------------------------- c c dimension r(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine fun (y,dy) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension y(1),dy(1),r(90) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c do 1 k=1,nr 1 r(k)= p(k) do 2 is=1,nisa 2 r(isa(1,is)) = r(isa(1,is))*y(isa(2,is))**isa(3,is) do 3 i=1,ny 3 dy(i) = 0. do 4 is=1,nisa 4 dy(isa(2,is))=dy(isa(2,is)) - dble(isa(3,is))*r(isa(1,is)) do 5 is=1,nisb 5 dy(isb(2,is))=dy(isb(2,is)) + sb(is)*r(isb(1,is)) return end c c *** KINAL ****************************************************** c * * c * JFY computes the Jacobian with respect to variables * c * * c **************************************************************** c c y - IN: actual concentrations c f - OUT: Jacobian matrix c fy(i,j) = d f(i) / d y(j) c ndy - IN: leading dimension of matrix fy c r - working area c c ---------------------------------------------------------------- c c dimension r(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine jfy(ndy,y,fy) implicit integer*2(i-n) implicit real*8(a-h,o-z) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) dimension fy(ndy,1),y(1),r(90) c do 1 i=1,ny do 1 j=1,ny 1 fy(i,j)=0. do 6 j=1,ny do 2 k=1,nr 2 r(k) = 0. do 3 is=1,nisa 3 if(isa(2,is).eq.j) r(isa(1,is)) = p(isa(1,is)) do 4 is=1,nisa if(r(isa(1,is)).eq.0.) goto 4 if(isa(2,is).eq.j) then r(isa(1,is)) = r(isa(1,is))*dble(isa(3,is)) 1 *y(isa(2,is))**(isa(3,is)-1) else r(isa(1,is)) = r(isa(1,is))*y(isa(2,is))**isa(3,is) endif 4 continue do 5 is=1,nisa if (r(isa(1,is)).eq.0.) goto 5 fy(isa(2,is),j)=fy(isa(2,is),j)-dble(isa(3,is))*r(isa(1,is)) 5 continue do 6 is=1,nisb if (r(isb(1,is)).eq.0.) goto 6 fy(isb(2,is),j)=fy(isb(2,is),j) + sb(is) * r(isb(1,is)) 6 continue return end c c *** KINAL ****************************************************** c * * c * JFP computes the Jacobian with respect to parameters * c * * c **************************************************************** c c y - IN: actual concentrations c fp - OUT: Jacobian matrix c fp(i,k) = d f(i) / d p(k) c ndy - IN: leading dimension of matrix fp c inp - IN: serial number of parameters , c related to which sensitivities are to be computed c r - working area c c ---------------------------------------------------------------- c c dimension r(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine jfp (np,ndy,y,fp,inp) implicit integer*2(i-n) implicit real*8(a-h,o-z) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) dimension y(1),fp(ndy,1),r(90),inp(1) c do 1 k=1,nr 1 r(k)= 1. do 2 is=1,nisa 2 r(isa(1,is))= r(isa(1,is)) * y(isa(2,is)) ** isa(3,is) do 3 i=1,ny do 3 k=1,np 3 fp(i,k)= 0. do 4 is=1,nisa do 4 k=1,np 4 if(inp(k).eq.isa(1,is)) 1 fp(isa(2,is),k)= - dble(isa(3,is)) * r(isa(1,is)) do 5 is=1,nisb do 5 k=1,np 5 if(inp(k).eq.isb(1,is)) 1 fp(isb(2,is),k)= sb(is) * r(isb(1,is)) return end c c *** KINAL ****************************************************** c * * c * QSA computes the quasi stationary sensitivity matrix * c * * c **************************************************************** c c dimension ip(ny) c subroutine qsa (ny,np,ndy,y,s,fy,plam,inp) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension s(ndy,1),fy(ny,1),inp(1),ip(50) c call jfp(np,ndy,y,s,inp) call jfy(ny,y,fy) do 1 i=1,ny fy(i,i)= fy(i,i) - plam 1 continue c 2 call dec(ny,fy,ip,ier) if(ier.ne.0) then write(6,10) ier stop endif c do 3 i=1,np call sol(ny,fy,ip,s(1,i)) 3 continue c 10 format(//5x,'DEC failed in QSA , error code: ',i3/ 1 //5x,'--- PROGRAM TERMINATED ---'///) return end c c *** KINAL ****************************************************** c * * c * SOL computes the solution of linear system A*X=B * c * * c **************************************************************** c c see C.B. Moler , Algorithm 423 , CACM 15 (1972) c subroutine sol(n,a,ip,b) integer*2 n,ip(1),nm1,k,kp1,m,i,kb,km1 double precision a(n,1),b(1),t if (n.eq.1) goto 50 nm1=n-1 do 20 k=1,nm1 kp1=k+1 m=ip(k) t=b(m) b(m)=b(k) b(k)=t do 10 i=kp1,n 10 b(i)=b(i)+a(i,k)*t 20 continue do 40 kb=1,nm1 km1=n-kb k=km1+1 b(k)=b(k)/a(k,k) t=-b(k) do 30 i=1,km1 30 b(i)=b(i)+a(i,k)*t 40 continue 50 b(1)=b(1)/a(1,1) return end c c *** KINAL ****************************************************** c * * c * DEC makes the Gaussian decomposition of matrix A * c * * c **************************************************************** c c see C.B. Moler , Algorithm 423 , CACM 15 (1972) c subroutine dec(n,a,ip,ier) integer*2 n,ip(1),ier,nm1,k,kp1,m,i,j double precision a(n,1),t ier=0 ip(n)=1 if(n.eq.1) goto 70 nm1=n-1 do 60 k=1,nm1 kp1=k+1 m=k do 10 i=kp1,n 10 if(dabs(a(i,k)).gt.dabs(a(m,k))) m=i ip(k)=m t=a(m,k) if(m.eq.k) goto 20 ip(n)=-ip(n) a(m,k)=a(k,k) a(k,k)=t 20 if (t.eq.0.d00) goto 80 t=1.d00/t do 30 i=kp1,n 30 a(i,k)=-a(i,k)*t do 50 j=kp1,n t=a(m,j) a(m,j)=a(k,j) a(k,j)=t if(t.eq.0.d00) goto 50 do 40 i=kp1,n 40 a(i,j)=a(i,j)+a(i,k)*t 50 continue 60 continue 70 k=n if (a(n,n).eq.0.d00) goto 80 return 80 ier=k ip(n)=0 return end c c *** KINAL ****************************************************** c * * c * RATE computes reaction rates * c * * c **************************************************************** c c y - IN: actual concentrations c r - OUT: reaction rates c c ---------------------------------------------------------------- c c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine rate(y,r) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension y(1),r(1) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c do 1 k=1,nr 1 r(k)= p(k) do 2 is=1,nisa r(isa(1,is)) = r(isa(1,is))*y(isa(2,is))**isa(3,is) 2 continue return end c c *** KINAL ****************************************************** c * * c * RTIME reads date and time * c * * c **************************************************************** c c This subroutine is machine and compiler dependent. c Recent version conforms to MS-FORTRAN 4.0 c subroutine rtime(cd,ct) implicit integer*2 (i-n) character*10 cd,ct c c date c call getdat(iy,im,id) iy1=iy/100 iy=iy-iy1*100 iy1=iy/10 iy2=iy-iy1*10 +48 iy1=iy1 +48 im1=im/10 im2=im-im1*10 +48 im1=im1 +48 id1=id/10 id2=id-id1*10 +48 id1=id1 +48 cd=char(im1)//char(im2)//'-'//char(id1)//char(id2)//'-' 1 //char(iy1)//char(iy2) c c time c call gettim(ih,im,is,ic) ih1=ih/10 ih2=ih-ih1*10 +48 ih1=ih1 +48 im1=im/10 im2=im-im1*10 +48 im1=im1 +48 is1=is/10 is2=is-is1*10 +48 is1=is1 +48 ct=char(ih1)//char(ih2)//':'//char(im1)//char(im2)//':' 1 //char(is1)//char(is2) return end