Re: OMEGA and SIGMA in INFN
From: Mats Karlsson <mats.karlsson@biof.uu.se>
Subject: Re: OMEGA and SIGMA in INFN
Date: 23 Jul 1997 17:16:23 -0400
Dear Mark,
Here is an INFN-routine I used quite some time ago to summarize the result from a 100 subproblem NONMEM run. Unfortunately, it is not commented and the coding is "ugly" but maybe of some use still.
Best regards,
Mats Karlsson
subroutine infn (icall,theta,datrec,indxs)
dimension theta(*),datrec(*),indxs(*)
dimension s(100),ss(100),sss(100)
dimension the1(100),the2(100),the3(100),om4(100)
dimension om1(100),om2(100),om3(100),si1(100)
dimension sr(100),ssr(100),e(100),sd(100)
double precision theta,object,s,ss,sss
double precision thetaf,omegaf,sigmaf
double precision the1,the2,the3
double precision om1,om2,om3,om4,si1
double precision sr,ssr
integer nprob,iprob,un,un2,un3
DATA un /80/
DATA un2 /90/
DATA un3 /95/
common /rocm6/ thetaf(20),omegaf(10,10),sigmaf(10,10)
common /rocm8/ object
common /rocm9/ iere,ierc
common /rocm10/ nrep,irep
if (icall.eq.0) then
nprob=100
ne=0
iprob=0
do 3 i=1,8
ssr(i)=0
3 sr(i)=0
do 5 i=1,nprob
s(i)=0
ss(i)=0
sss(i)=0
the1(i)=0
the2(i)=0
the3(i)=0
om4(i)=0
om1(i)=0
om2(i)=0
om3(i)=0
5 si1(i)=0
endif
if (icall.eq.1) then
iprob=iprob+1
endif
if (icall.eq.3) then
s(iprob)=object
ss(iprob)=iere
sss(iprob)=iprob
the1(iprob)=thetaf(1)
the2(iprob)=thetaf(2)
the3(iprob)=thetaf(3)
om1(iprob)=omegaf(1,1)
om2(iprob)=omegaf(2,2)
om3(iprob)=omegaf(3,3)
om4(iprob)=omegaf(4,4)
si1(iprob)=sigmaf(1,1)
if(iprob.eq.nprob) then
do 6 i=1,nprob
6 write(un,*)s(i),ss(i),sss(i)
if (iprob.eq.nprob) then
do 7 i=1,nprob
7 write(un2,8)the1(i),the2(i),the3(i),
1om1(i),om2(i),om3(i),om4(i),si1(i),sss(i)
8 FORMAT(9F9.4X)
endif
endif
if (iere.eq.0) then
ne=ne+1
sr(1)=sr(1)+the1(iprob)
ssr(1)=ssr(1)+the1(iprob)**2
sr(2)=sr(2)+the2(iprob)
ssr(2)=ssr(2)+the2(iprob)**2
sr(3)=sr(3)+the3(iprob)
ssr(3)=ssr(3)+the3(iprob)**2
sr(4)=sr(4)+om4(iprob)
ssr(4)=ssr(4)+om4(iprob)**2
sr(5)=sr(5)+om1(iprob)
ssr(5)=ssr(5)+om1(iprob)**2
sr(6)=sr(6)+om2(iprob)
ssr(6)=ssr(6)+om2(iprob)**2
sr(7)=sr(7)+om3(iprob)
ssr(7)=ssr(7)+om3(iprob)**2
sr(8)=sr(8)+si1(iprob)
ssr(8)=ssr(8)+si1(iprob)**2
endif
if (iprob.eq.nprob) then
if (ne.gt.0) then
e(1)=sr(1)/ne
e(2)=sr(2)/ne
e(3)=sr(3)/ne
e(4)=sr(4)/ne
e(5)=sr(5)/ne
e(6)=sr(6)/ne
e(7)=sr(7)/ne
e(8)=sr(8)/ne
if (ne.gt.1) then
sd(1)=sqrt((ssr(1)-sr(1)**2/ne)/(ne-1))
sd(2)=sqrt((ssr(2)-sr(2)**2/ne)/(ne-1))
sd(3)=sqrt((ssr(3)-sr(3)**2/ne)/(ne-1))
sd(4)=sqrt((ssr(4)-sr(4)**2/ne)/(ne-1))
sd(5)=sqrt((ssr(5)-sr(5)**2/ne)/(ne-1))
sd(6)=sqrt((ssr(6)-sr(6)**2/ne)/(ne-1))
sd(7)=sqrt((ssr(7)-sr(7)**2/ne)/(ne-1))
sd(8)=sqrt((ssr(8)-sr(8)**2/ne)/(ne-1))
endif
endif
do 9 i=1,8
9 write(un3,*)ne,e(i),sd(i)
endif
endif
end