program bootpi implicit none integer i,j,k,nsubj,nrespin,nboot,nmax,count,nfourtrm,kr,seed integer addpts,totaddpts,rowsA integer maxnrespin,maxsubj,maxfourtrm,maxboot,maxtotaddpts parameter (maxnrespin = 200,maxsubj = 70,maxfourtrm = 50, + maxboot = 1000,maxtotaddpts = 55) integer maxrow,maxnmax parameter (maxrow = maxnrespin + maxtotaddpts, + maxnmax = maxboot * maxsubj) C ************************************************* C Name: nspi C Purpose: main program C Date Updated: 2-22-99 C Change Record: C When Who What C Previous 2-22 BW created C 2-22/2-23 JJP implemented batchfile processing C 3-31 JJP began multiple job logic C 11/27/01 LoriLyn Price & TJS C add'l diagnostics added C ************************************************* c c METHOD c c This program fits a Fourier model to the data for each subject. c A mean curve is estimated by taking the sample average of c the fitted values over all subjects at each input point. The c sample covariance matrix of the parameter estimates is also c obtained. From this matrix, an estimate of variability is c calculated at each input point. A bootstrap procedure is c then applied to estimate the critical value needed for the c construction of simultaneous prediction bands. These bands c are calculated using the aforementioned estimates of the c mean, variance and critical value. c c DEFINITIONS c c maxnrespin = MAXIMUM number of points per curve on input c maxsubj = MAXIMUM number of curves to be analyzed c maxfourtrm = MAXIMUM number of Fourier terms c maxboot = MAXIMUM number of bootstrap reps c maxtotaddpts = MAXIMUM number of points for curve extension c (total for both sides) c c Note: The numerical values of these parameters can be c modified for use of this program on a PC. c c INPUT c c The program will prompt the user for the name of the parameter c and data input files. c c The parameter input file should contain numerical values for c the following parameters in column format: c nsubj = number of curves to be analyzed c nrespin = number of points per curve on input c nfourtrm = number of Fourier terms c addpts = number of points for each side of curve c extension c nboot = number of bootstrap reps c covprob = coverage probability for prediction intervals c seed = non-zero integer used for random number generation c (seed must be between -32765 and 32765 excluding 0) c c The data input file should contain the curves to be c analyzed in column format (i.e. columns represent subjects). c c OUTPUT c c The program will prompt the user for the name of the output c file. c c The output file will contain four columns. The first c column gives the points at which the mean curve and the c prediction limits were evaluated. The last three columns c give the values of the lower prediciton limit, mean curve, c and upper prediction limit respectively. c c The output file will also contain the maximum c absolute deviation between the observed and fitted c response for each subject. The maximum overall c absolute deviation is also provided. integer jpvt(maxfourtrm) real*8 tol,m,covprob,quantile real*8 maxdiff real*8 favg(maxnrespin),favgbs(maxnrespin),favgls(maxfourtrm), + favglsb(maxfourtrm) real*8 theta(maxrow),qraux(maxfourtrm),work(maxfourtrm) real*8 sigma(maxnrespin),sigmabs(maxnrespin) real*8 RSD(maxrow),B(maxrow),X(maxfourtrm) real*8 maxdev(maxnmax),predu(maxnrespin),predl(maxnrespin) real*8 augtheta(maxtotaddpts),deriv(maxsubj) real*8 diff(maxsubj) real*8 temp(maxnrespin),difftemp(maxrow) integer sam(maxboot,maxsubj) real*8 A(maxrow,maxfourtrm),Y(maxrow,maxsubj), + ls(maxfourtrm,maxsubj),f(maxnrespin,maxsubj), + covhat(maxfourtrm,maxfourtrm),bf(maxnrespin,maxsubj), + lsbs(maxfourtrm,maxsubj),covhatbs(maxfourtrm,maxfourtrm), + dev(maxnrespin,maxsubj) real*8 Atemp(maxrow,maxfourtrm),Ytemp(maxrow,maxsubj) real*8 lsc(maxfourtrm,maxsubj),lsbsc(maxfourtrm,maxsubj) real*8 mat(maxnrespin,maxfourtrm) real*8 currmin, currmax, range, relrange Character*40 doutfile,coutfile,crunname,cBatchFile,datfile Integer NumArgs, iHold,linecounter, seed_hold,iargc Logical lAbort,lFinished lFinished = .false. NumArgs = iargc() iHold = 1 if(NumArgs.ne.1)then Write(6,*)"Name of Batch file was not entered." Write(6,*)"Please enter name of Batch File: " Read(5,*)cBatchFile else Call Getarg(iHold,cBatchFile) end if linecounter = 0 OPEN(UNIT=20,FILE=cBatchFile,status="old") do while(lFinished .eq. .false.) 1 call readdata(nsubj,nrespin,nfourtrm,addpts,nboot,covprob,seed, + Y,maxsubj,maxrow,doutfile,coutfile,cRunName,lAbort, * lFinished,linecounter,datfile) if(lAbort.eq..true.)then go to 1 end if if(lFinished.eq..true.)then go to 999 end if seed_hold = seed totaddpts = 2*addpts rowsA = nrespin + 2*addpts currmin = Y(1,1) currmax = Y(1,1) do 55 i=1,nsubj do 56 j=1,nrespin currmin=min(currmin, Y(j,i)) currmax=max(currmax, Y(j,i)) 56 continue 55 continue call calctheta(theta,nrespin) call calcaugtheta(augtheta,addpts,theta,nrespin,totaddpts) call derivative(nrespin,nsubj,Y,theta,deriv,maxrow,1,Ytemp) call polyint(nsubj,nrespin,addpts,deriv,theta,augtheta,Y, + maxrow,totaddpts,1,Ytemp) call derivative(nrespin,nsubj,Y,theta,deriv,maxrow,2,Ytemp) call polyint(nsubj,nrespin,addpts,deriv,theta,augtheta,Y, + maxrow,totaddpts,2,Ytemp) call des(theta,nrespin,nfourtrm,A,maxrow) call augdes(augtheta,addpts,nrespin,nfourtrm,A,maxrow, * totaddpts,Atemp) tol=1.d-15 call DQRANK(A,maxrow,rowsA,nfourtrm,tol,kr,jpvt,qraux,work) call lsest(A,maxrow,rowsA,nfourtrm,nsubj,kr,B,X,RSD,jpvt,qraux, + Y,ls,maxfourtrm) call calctheta(theta,nrespin) call des(theta,nrespin,nfourtrm,A,maxrow) call func(nrespin,nfourtrm,nsubj,ls,maxfourtrm,f,maxnrespin,A, + maxrow) call diagnostic(nrespin,nsubj,Ytemp,f,diff,difftemp,maxrow, + maxnrespin) call fbar(nrespin,nsubj,favg,1,ls,lsbs,maxfourtrm,f,bf, * maxnrespin) call fbar(nfourtrm,nsubj,favgls,3,ls,lsbs,maxfourtrm,f,bf, + maxnrespin) call cov(nfourtrm,nsubj,favgls,3,ls,covhat,lsbs,covhatbs, + maxfourtrm,lsc,lsbsc) call sig(sigma,nrespin,nfourtrm,5,covhat,covhatbs,maxfourtrm,A, + maxrow,mat,maxnrespin) nmax = nboot * nsubj count = 0 call bsample(nboot,nsubj,seed,sam,maxboot) do 10 i=1,maxnmax maxdev(i)=-1.d0 10 continue do 20 k=1,nboot call bstrapf(nrespin,nsubj,k,f,bf,maxnrespin,sam,maxboot) call fbar(nrespin,nsubj,favgbs,2,ls,lsbs,maxfourtrm, + f,bf,maxnrespin) call lsestbs(nfourtrm,nsubj,k,ls,lsbs,maxfourtrm,sam, * maxboot) call fbar(nfourtrm,nsubj,favglsb,4,ls,lsbs,maxfourtrm, + f,bf,maxnrespin) call cov(nfourtrm,nsubj,favglsb,4,ls,covhat,lsbs, + covhatbs,maxfourtrm,lsc,lsbsc) call sig(sigmabs,nrespin,nfourtrm,6,covhat,covhatbs, * maxfourtrm,A,maxrow,mat,maxnrespin) call norm(nrespin,nsubj,favgbs,sigmabs,f,dev,maxnrespin) call maximize(maxdev,nrespin,nsubj,count,nmax,dev, + maxnrespin,temp) count = count + nsubj 20 continue call sort(nmax,maxdev) m = quantile(maxdev,nmax,covprob) do 30 i=1,nrespin predu(i) = favg(i) + m * sigma(i) predl(i) = favg(i) - m * sigma(i) 30 continue open(15, file=coutfile,status='unknown') open(16, file=doutfile,status='unknown') write(16,*)"Input File: ",datfile write(16,*)"Confidence Limit Output File: ",coutfile Write(16,*)"Seed: ",seed_hold write(16,*) " " write(16,210) write(16,211) covprob * 100.d0 write(16,212) nboot write(16,*) 210 format(1x,'Simultaneous Prediction Intervals') 211 format(1x,'Coverage Probability is ',f16.4,' %') 212 format(1x,'Number of Bootstrap Reps = ',I10) do 40 i=1,nrespin write(15,213) theta(i),predl(i),favg(i),predu(i) 40 continue 213 format(1x,4(f16.5,1x)) write(16,214) nfourtrm write(16,215) addpts 214 format(1x,'Number of Fourier Terms = ',I10) 215 format(1x,'Number of Curve Extension Points = ',I10) write(16,*) write(16,216) write(16,217) write(16,*) do 50 i=1,nsubj write(16,218) i, diff(i) 50 continue call sort(nsubj,diff) maxdiff = diff(nsubj) range = currmax - currmin relrange = maxdiff/range write(16,*) write(16,219) maxdiff write(16,220) relrange 216 format(1x,'Maximum absolute deviation between observed and') 217 format(1x,'fitted response for each curve:') 218 format(1x,i5,3x,f16.5) 219 format(1x,'Maximum absolute deviation over curves = ',f16.5) 220 format(1x,'Maximum relative deviation over curves = ',f16.5) close(15) close(16) write(6,*)" " end do 999 Close(20) end subroutine readdata(nsubj,nrespin,nfourtrm,addpts,nboot,covprob, + seed,Y,maxsubj,maxrow,doutfile,coutfile,cRunName,lAbort, * lFinished,linecounter,datfile) C ************************************************* C Name: readdata C Purpose: read in the data from the batchfile C Date Updated: 2-22-99 C Change Record: C When Who What C Previous 2-22 BW created C 2-22/2-23 JJP implemented batchfile processing C 4-1 JJP began other input method batch file options C ************************************************* implicit none integer i,j,nsubj,nrespin,nfourtrm,addpts,nboot,seed integer maxsubj,maxrow,linecounter,iInputMethod,iNumLinesOmit integer iColumnNumber,internalcount,location1 integer location2,iPoint,iSubject,jk character*40 datfile,doutfile,coutfile,cRunName,cSubDirFile character*64 cInString,cinStringReg,Tranlc,cDirectoryName character*64 cJunkLine,cOpenFile real*8 covprob real*8 Y(maxrow,maxsubj),XCoord,Dimension(3) Logical lDataFile, lDOutFile, lNumCurves, lNumPoints, lNumTerms Logical lNumPointsExt, lNumBoot, lConv, lSeed, lFirst, lAbort Logical lDone, lCOutFile, lFinished,lSubDirFile Logical lRun,lPrint cRunName = '' cSubDirFile = '' datfile = '' doutfile = '' coutfile = '' iInputMethod = 1 iNumLinesOmit = 0 iColumnNumber = 0 internalcount = 0 lDataFile = .false. lDOutFile = .false. lCOutfile = .false. lNumCurves = .false. lNumPoints = .false. lNumTerms = .false. lNumPointsExt = .false. lNumBoot = .false. lSubDirFile = .false. lConv =.false. lSeed = .false. lDone = .false. lAbort = .false. lRun = .false. lPrint = .false. do while (lDone.eq..false.) internalcount = internalcount + 1 if((internalcount.gt.1).and.(lPrint.eq..false.).and. * (lRun.eq..false.))then write(6,*)"**-> Begin Job" lPrint = .true. end if linecounter = linecounter + 1 cInString='' cInStringReg='' Read(20,'(A64)') cInString Call Squeeze(cInString) cInStringReg=Tranlc(cInString) if(cInString(1:8).eq.'runname=')then cRunName = cInStringReg(9:) lRun = .true. write(6,*)"**-> Begin ",cRunName else if(cInString(1:14).eq.'datainputfile=')then datfile = cInStringReg(15:) lDataFile = .true. else if(cInString(1:21).eq.'diagnosticoutputfile=')then doutfile = cInStringReg(22:) lDOutFile = .true. else if(cInString(1:20).eq.'conflimitoutputfile=')then coutfile = cInStringReg(21:) lCOutFile = .true. else if(cInString(1:10).eq.'numcurves=')then Read(cInStringReg(11:),*)nsubj lNumCurves = .true. if(nsubj.lt.1)then lNumCurves = .false. end if else if(cInString(1:18).eq.'numpointspercurve=')then Read(cInStringReg(19:),*)nrespin lNumPoints = .true. if(nrespin.lt.1)then lNumPoints = .false. end if else if(cInString(1:16).eq.'numfourierterms=')then Read(cInStringReg(17:),*)nfourtrm lNumTerms = .true. if(nfourtrm.lt.1)then lNumTerms = .false. end if else if(cInString(1:19).eq.'numpointsextension=')then Read(cInStringReg(20:),*)addpts lNumPointsExt = .true. if(addpts.lt.1)then lNumPointsExt = .false. end if else if(cInString(1:17).eq.'numbootstrapreps=')then Read(cInStringReg(18:),*)nboot lNumBoot = .true. if(nboot.lt.1)then lNumBoot = .false. end if else if(cInString(1:8).eq.'covprob=')then Read(cInStringReg(9:),*)covprob lConv = .true. if((covprob.le.0).or.(covprob.ge.1))then lConv = .false. end if else if(cInString(1:5).eq.'seed=')then Read(cInStringReg(6:),*)seed lSeed = .true. if((seed.lt.-32765).or.(seed.gt.32765).or. * (seed.eq.0))then lSeed = .false. end if else if(cInString(1:16).eq.'datainputmethod=')then Read(cInStringReg(17:),*)iInputMethod else if(cInString(1:15).eq.'omitnumoflines=')then if(iInputMethod.ne.2)then write(6,*)"OmitNumberOfLines was read from the batch "// * "file," write(6,*)" but DataInputMethod was not set to 2." write(6,*)"Therefore, OmitNumberOfLines will not be "// * "used." write(6,*)"--" end if Read(cInStringReg(16:),*)iNumLinesOmit else if(cInString(1:16).eq.'columntoanalyze=')then if(iInputMethod.ne.2)then write(6,*)"ColumnToAnalyze was read from the batch "// * "file," write(6,*)" but DataInputMethod was not set to 2." write(6,*)"Therefore, ColumnToAnalyze will not be "// * "used." write(6,*)"--" end if Read(cInStringReg(17:),*)iColumnNumber else if(cInString(1:17).eq.'subdirectoryfile=')then if(iInputMethod.ne.2)then write(6,*)"SubDirectoryFile was read from the batch "// * "file," write(6,*)" but DataInputMethod was not set to 2." write(6,*)"Therefore, SubDirectoryFile will not be "// * "used." write(6,*)"--" end if cSubDirFile = cInStringReg(18:) lSubDirFile = .true. else if(cInString(1:4).eq.'stop')then lDone = .true. else if(cInString(1:3).eq.'end')then lFinished = .true. go to 21 else if((cInString(1:1).eq.' ').or. * (cInString(1:1).eq.'#'))then cInString = '' internalcount = internalcount - 1 else Write(6,*)'Unknown error in batch file at '// * 'line ', lineCounter, '.' Write(6,*)'Error ignored. Check spelling '// * 'and input format.' write(6,*)"--" end if end do lFirst = .true. if (lDataFile.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"DataInputFile" end if if (lDOutFile.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"DiagnosticOutputFile" end if if (lCOutFile.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"ConfLimitOutputFile" end if if (lNumCurves.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"NumCurves" end if if (lNumPoints.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"NumPointsPerCurve" end if if (lNumTerms.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"NumFourierTerms" end if if (lNumPointsExt.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"NumPointsExtension" end if if (lNumBoot.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"NumBootstrapReps" end if if (lConv.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"CovProb" end if if (lSeed.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"Seed" end if if(iInputMethod.eq.2)then if((iColumnNumber.lt.2).or.(iColumnNumber.gt.4))then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the "// * "batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"ColumnToAnalyze" end if if(lSubDirFile.eq..false.)then if(lFirst.eq..true.)then Write(6,*)"The following parameters are missing" Write(6,*)" or the values are incorrect in the "// * "batchfile." Write(6,*)"Check the name and spelling of the "// * "parameter names." Write(6,*)"--------------------" lFirst = .false. end if Write(6,*)"SubDirectoryFile" end if end if 19 if(lFirst.eq..false.)then write(6,*)" " write(6,*)"This job will now abort. Fix batch file "// * "and rerun." Write(6,*)" " lAbort = .true. Return end if C C Read in input file(s) C if(iInputMethod.eq.1)then open(17,file=datfile,status='old') rewind 17 do 20 i=1,nrespin read(17,*) (Y(i,j),j=1,nsubj) 20 continue close(17) return else if(iInputMethod.eq.2)then Open(18,file=cSubDirFile,status='old') location2 = index(datfile,' ') do iSubject=1,nsubj Read(18,'(A64)')cDirectoryName location1 = index(cDirectoryName,' ') if(cDirectoryName(location1-1:location1-1).eq.'/')then location1 = location1 - 1 end if cOpenFile = cDirectoryName(1:location1-1) // '/' // * datfile(1:location2-1) Open(19,file=cOpenFile,status='old') if (iNumLinesOmit.eq.0)go to 201 do jk=1,iNumLinesOmit Read(19,'(a64)')cJunkLine end do 201 Continue do iPoint=1,nrespin Read(19,*)XCoord,(Dimension(i),i=1,3) Y(iPoint,iSubject) = Dimension(iColumnNumber) end do Close(19) end do else if(iInputMethod.gt.2)then Write(6,*)"DataInputMethod can only be a 1 or 2." lFirst = .false. go to 19 end if 21 return end subroutine calctheta(theta,dim) implicit none integer i,dim real*8 theta(dim) do 10 i=1,dim theta(i)=dfloat(100*(i-1))/dfloat(3*dim) 10 continue return end subroutine calcaugtheta(augtheta,addpts,theta,nrespin,totaddpts) implicit none integer i,addpts,nrespin,totaddpts real*8 theta1 real*8 augtheta(totaddpts),theta(nrespin) do 10 i=1,addpts theta1 = theta(nrespin) augtheta(addpts-i+1)=-theta1 * (dfloat(i)/dfloat(addpts)) augtheta(addpts+i)=theta1 * + (1.d0 + dfloat(i)/dfloat(addpts)) 10 continue return end subroutine derivative(nrespin,nsubj,Y,theta,deriv,maxrow,eind, + Ytemp) implicit none integer i,j,nrespin,nsubj,eind integer maxrow real*8 A1,A2,A3,B1,B2 real*8 x(3),z(3) real*8 theta(nrespin),deriv(nsubj) real*8 Y(maxrow,nsubj),Ytemp(maxrow,nsubj) if (eind.eq.1) then do 10 i=1,nsubj do 20 j=1,3 x(j)=theta(j) z(j)=Y(j,i) 20 continue A1=z(1)/((x(1)-x(2))*(x(1)-x(3))) A2=z(2)/((x(2)-x(1))*(x(2)-x(3))) A3=z(3)/((x(3)-x(1))*(x(3)-x(2))) B2=A1 * (x(2)+x(3)) + A2 * (x(1)+x(3)) + A3 * (x(1)+x(2)) deriv(i)=-B2 10 continue elseif (eind.eq.2) then do 30 i=1,nsubj do 40 j=1,3 x(4-j)=theta(nrespin-(j-1)) z(4-j)=Ytemp(nrespin-(j-1),i) 40 continue A1=z(1)/((x(1)-x(2))*(x(1)-x(3))) A2=z(2)/((x(2)-x(1))*(x(2)-x(3))) A3=z(3)/((x(3)-x(1))*(x(3)-x(2))) B1=A1 + A2 + A3 B2=A1 * (x(2)+x(3)) + A2 * (x(1)+x(3)) + A3 * (x(1)+x(2)) deriv(i)=2.d0 * B1 * theta(nrespin) - B2 30 continue endif return end subroutine polyint(nsubj,nrespin,addpts,deriv,theta,augtheta, + Y,maxrow,totaddpts,eind,Ytemp) implicit none integer i,j,nsubj,nrespin,addpts integer maxrow,totaddpts,eind real*8 u1,u2,u3,p0,p1,p2,theta1,theta2 real*8 deriv(nsubj),theta(nrespin),augtheta(totaddpts) real*8 Y(maxrow,nsubj),Ytemp(maxrow,nsubj) if (eind.eq.1) then do 10 j=1,nsubj do 20 i=1,nrespin Ytemp(i,j)=Y(i,j) 20 continue 10 continue do 30 j=1,nsubj do 40 i=1,nrespin Y(addpts+i,j)=Ytemp(i,j) 40 continue 30 continue theta1 = theta(nrespin) do 50 j=1,nsubj p0 = Ytemp(1,j) p1 = deriv(j) p2 = p1/theta1 do 60 i=1,addpts theta2 = augtheta(i) Y(i,j) = p0 + p1 * theta2 + p2 * (theta2)**2.d0 60 continue 50 continue elseif (eind.eq.2) then theta1 = theta(nrespin) do 70 j=1,nsubj u1 = deriv(j) u2 = Ytemp(1,j) u3 = Ytemp(nrespin,j) p0 = -2.d0 * theta1 * u1 + u2 p1 = 3.d0 * u1 + 2.d0 * (u3 - u2)/theta1 p2 = -u1/theta1 + (u2 - u3)/((theta1)**2.d0) do 80 i=1,addpts theta2 = augtheta(addpts+i) Y(addpts+nrespin+i,j) = p0 + p1 * theta2 + + p2 * (theta2)**2.d0 80 continue 70 continue endif return end subroutine des(theta,dim,nfourtrm,A,maxrow) implicit none integer i,j,c,nfourtrm,dim integer maxrow real*8 pi,k,theta(dim),d1,d2 real*8 A(maxrow,nfourtrm) data pi /3.141592654d0/ k = pi/50.d0 do 10 j=1,nfourtrm c=mod(j,2) do 20 i=1,dim if (j.eq.1) then A(i,j) = 1.d0 elseif (c.eq.0) then d1 = k*dfloat(j/2)*theta(i) A(i,j) = dcos(d1) elseif (j.gt.1 .and. c.eq.1) then d2 = k*dfloat((j-1)/2)*theta(i) A(i,j) = dsin(d2) endif 20 continue 10 continue return end subroutine augdes(augtheta,addpts,nrespin,nfourtrm,A,maxrow, + totaddpts,Atemp) implicit none integer i,j,c,addpts,nrespin,nfourtrm integer maxrow,totaddpts real*8 pi,k,augtheta(totaddpts),d1,d2 real*8 A(maxrow,nfourtrm),Atemp(maxrow,nfourtrm) data pi /3.141592654d0/ k = pi/50.d0 do 10 j=1,nfourtrm do 20 i=1,nrespin Atemp(i,j) = A(i,j) 20 continue 10 continue do 30 j=1,nfourtrm do 40 i=1,nrespin A(addpts+i,j) = Atemp(i,j) 40 continue 30 continue do 50 j=1,nfourtrm c=mod(j,2) do 60 i=1,addpts if (j.eq.1) then A(i,j) = 1.d0 elseif (c.eq.0) then d1 = k*dfloat(j/2)*augtheta(i) A(i,j) = dcos(d1) elseif (j.gt.1 .and. c.eq.1) then d2 = k*dfloat((j-1)/2)*augtheta(i) A(i,j) = dsin(d2) endif 60 continue 50 continue do 70 j=1,nfourtrm c=mod(j,2) do 80 i=1,addpts if (j.eq.1) then A(addpts+nrespin+i,j) = 1.d0 elseif (c.eq.0) then d1 = k*dfloat(j/2)*augtheta(addpts+i) A(addpts+nrespin+i,j) = dcos(d1) elseif (j.gt.1 .and. c.eq.1) then d2 = k*dfloat((j-1)/2)*augtheta(addpts+i) A(addpts+nrespin+i,j) = dsin(d2) endif 80 continue 70 continue return end subroutine lsest(A,maxrow,rowsA,nfourtrm,nsubj,kr,B,X,RSD, + jpvt,qraux,Y,ls,maxfourtrm) implicit none integer i,j,k,rowsA,nfourtrm,kr,nsubj,jpvt(nfourtrm) integer maxfourtrm,maxrow real*8 qraux(nfourtrm) real*8 A(maxrow,nfourtrm),RSD(rowsA),B(rowsA),X(nfourtrm) real*8 Y(maxrow,nsubj),ls(maxfourtrm,nsubj) do 10 j=1,nsubj do 20 i=1,rowsA B(i) = Y(i,j) 20 continue call DQRLSS(A,maxrow,rowsA,nfourtrm,kr,B,X,RSD,jpvt,qraux) do 30 k=1,nfourtrm ls(k,j) = X(k) 30 continue 10 continue return end subroutine lsestbs(nfourtrm,nsubj,k,ls,lsbs,maxfourtrm,sam, + maxboot) implicit none integer i,j,k,nfourtrm,nsubj integer maxfourtrm,maxboot integer sam(maxboot,nsubj) real*8 ls(maxfourtrm,nsubj),lsbs(maxfourtrm,nsubj) do 10 j=1,nsubj do 20 i=1,nfourtrm lsbs(i,j)=ls(i,sam(k,j)) 20 continue 10 continue return end subroutine func(nrespin,nfourtrm,nsubj,ls,maxfourtrm,f,maxnrespin, + A,maxrow) implicit none integer i,j,k,nrespin,nfourtrm,nsubj integer maxfourtrm,maxnrespin,maxrow real*8 A(maxrow,nfourtrm),ls(maxfourtrm,nsubj),f(maxnrespin,nsubj) do 10 i=1,nrespin do 20 j=1,nsubj f(i,j) = 0.d0 20 continue 10 continue do 30 i=1,nrespin do 40 j=1,nsubj do 50 k=1,nfourtrm f(i,j)=f(i,j) + A(i,k) * ls(k,j) 50 continue 40 continue 30 continue return end subroutine bstrapf(nrespin,nsubj,k,f,bf,maxnrespin,sam,maxboot) implicit none integer i,j,k,nsubj,nrespin integer maxnrespin,maxboot integer sam(maxboot,nsubj) real*8 f(maxnrespin,nsubj),bf(maxnrespin,nsubj) do 10 j=1,nsubj do 20 i=1,nrespin bf(i,j)=f(i,sam(k,j)) 20 continue 10 continue return end subroutine fbar(nrow,nsubj,avg,ind,ls,lsbs,maxfourtrm, + f,bf,maxnrespin) implicit none integer i,j,nrow,nsubj,ind integer maxfourtrm,maxnrespin real*8 avg(nrow) real*8 ls(maxfourtrm,nsubj),f(maxnrespin,nsubj), + bf(maxnrespin,nsubj),lsbs(maxfourtrm,nsubj) do 10 i=1,nrow avg(i)=0.d0 10 continue do 20 i=1,nrow do 30 j=1,nsubj if (ind.eq.1) then avg(i)=avg(i) + f(i,j) elseif (ind.eq.2) then avg(i)=avg(i) + bf(i,j) elseif (ind.eq.3) then avg(i)=avg(i) + ls(i,j) elseif (ind.eq.4) then avg(i)=avg(i) + lsbs(i,j) endif 30 continue avg(i)=avg(i)/dfloat(nsubj) 20 continue return end subroutine cov(nfourtrm,nsubj,avg,ind,ls,covhat,lsbs, + covhatbs,maxfourtrm,lsc,lsbsc) implicit none integer i,j,k,nfourtrm,nsubj,ind integer maxfourtrm real*8 avg(nfourtrm),lsc(maxfourtrm,nsubj),lsbsc(maxfourtrm,nsubj) real*8 ls(maxfourtrm,nsubj),covhat(maxfourtrm,nfourtrm), + lsbs(maxfourtrm,nsubj),covhatbs(maxfourtrm,nfourtrm) do 10 i=1,nfourtrm do 20 j=1,nfourtrm if (ind.eq.3) then covhat(i,j)=0.d0 elseif (ind.eq.4) then covhatbs(i,j)=0.d0 endif 20 continue 10 continue do 30 i=1,nfourtrm do 40 j=1,nsubj if (ind.eq.3) then lsc(i,j)=ls(i,j) - avg(i) elseif (ind.eq.4) then lsbsc(i,j)=lsbs(i,j) - avg(i) endif 40 continue 30 continue do 50 i=1,nfourtrm do 60 j=1,i do 70 k=1,nsubj if (ind.eq.3) then covhat(i,j)=lsc(i,k) * lsc(j,k) + covhat(i,j) elseif (ind.eq.4) then covhatbs(i,j)=lsbsc(i,k) * lsbsc(j,k) + + covhatbs(i,j) endif 70 continue if (ind.eq.3) then covhat(i,j)=covhat(i,j)/dfloat(nsubj - 1) elseif (ind.eq.4) then covhatbs(i,j)=covhatbs(i,j)/dfloat(nsubj - 1) endif 60 continue 50 continue do 80 i=1,nfourtrm - 1 do 90 j=i+1,nfourtrm if (ind.eq.3) then covhat(i,j) = covhat(j,i) elseif(ind.eq.4) then covhatbs(i,j)=covhatbs(j,i) endif 90 continue 80 continue return end subroutine sig(stdev,nrespin,nfourtrm,ind,covhat,covhatbs, + maxfourtrm,A,maxrow,mat,maxnrespin) implicit none integer i,j,k,nrespin,nfourtrm,ind integer maxfourtrm,maxrow,maxnrespin real*8 mat(maxnrespin,nfourtrm),stdev(nrespin) real*8 A(maxrow,nfourtrm),covhat(maxfourtrm,nfourtrm), + covhatbs(maxfourtrm,nfourtrm) do 10 i=1,nrespin stdev(i)=0.d0 do 20 j=1,nfourtrm mat(i,j)=0.d0 20 continue 10 continue do 30 i=1,nrespin do 40 j=1,nfourtrm do 50 k=1,nfourtrm if(ind.eq.5) then mat(i,j)=A(i,k) * covhat(k,j) + mat(i,j) elseif (ind.eq.6) then mat(i,j)=A(i,k) * covhatbs(k,j) + mat(i,j) endif 50 continue 40 continue 30 continue do 60 i=1,nrespin do 70 j=1,nfourtrm stdev(i)=stdev(i) + mat(i,j) * A(i,j) 70 continue stdev(i)=dsqrt(stdev(i)) 60 continue return end subroutine norm(nrespin,nsubj,favgbs,sigmabs,f,dev,maxnrespin) implicit none integer i,j,nrespin,nsubj integer maxnrespin real*8 favgbs(nrespin),sigmabs(nrespin) real*8 f(maxnrespin,nsubj),dev(maxnrespin,nsubj) do 10 i=1,nrespin do 20 j=1,nsubj dev(i,j)=dabs((f(i,j) - favgbs(i))/sigmabs(i)) 20 continue 10 continue return end subroutine maximize(maxdev,nrespin,nsubj,count,nmax,dev, + maxnrespin,temp) implicit none integer i,j,nrespin,nsubj,count,nmax integer maxnrespin real*8 temp(nrespin),maxdev(nmax) real*8 dev(maxnrespin,nsubj) do 10 i=1,nsubj do 20 j=1,nrespin temp(j)=dev(j,i) 20 continue call sort(nrespin,temp) maxdev(i+count)=temp(nrespin) 10 continue return end subroutine bsample(nboot,nsubj,seed,sam,maxboot) implicit none integer i,j,p,seed,nboot,nsubj integer maxboot integer sam(maxboot,nsubj) real u,uni do 10 i=1,nboot do 20 j=1,nsubj u=uni(seed) seed=0 p=int(float(nsubj) * u) if (p.le.(nsubj - 1)) then sam(i,j)=p + 1 elseif (p.eq.nsubj) then sam(i,j)=nsubj endif 20 continue 10 continue return end function quantile(maxdev,nmax,covprob) implicit none integer i,nmax,q real*8 covprob real*8 j,delta,ratio,quantile,maxdev(nmax) do 10 i=1,nmax j=dfloat(i)/dfloat(nmax) if (j.eq.covprob) then quantile=maxdev(i) goto 200 elseif (j.gt.covprob .and. i.ge.2) then q=i-1 goto 100 elseif (j.gt.covprob .and. i.eq.1) then print *, 'Your coverage probability is too small!' print *, 'Try again!' stop endif 10 continue 100 delta=maxdev(q+1) - maxdev(q) ratio=dfloat(nmax) * covprob - dfloat(q) quantile=maxdev(q) + delta * ratio 200 return end subroutine diagnostic(nrespin,nsubj,Ytemp,f,diff,difftemp, + maxrow,maxnrespin) implicit none integer i,j,nrespin,nsubj integer maxrow,maxnrespin real*8 delta real*8 diff(nsubj) real*8 Ytemp(maxrow,nsubj),f(maxnrespin,nsubj) real*8 difftemp(nrespin) do 10 j=1,nsubj do 20 i=1,nrespin delta=Ytemp(i,j) - f(i,j) difftemp(i)=dabs(delta) 20 continue call sort(nrespin,difftemp) diff(j)=difftemp(nrespin) 10 continue return end c REFERENCE: Press, William H., Flannery, Brian P., Teukolsky, c Saul A. and Vetterling, William T., 1986, c Numerical Recipies: The Art of Scientific c Computing (Cambridge University Press). subroutine sort(n,ra) implicit none integer n,l,ir,i,j real*8 rra, ra(n) l=n/2 + 1 ir=n 10 continue if(l.gt.1) then l=l-1 rra=ra(l) else rra=ra(ir) ra(ir)=ra(1) ir=ir - 1 if(ir.eq.1) then ra(1)=rra return endif endif i=l j=l + l 20 if(j.le.ir) then if(j.lt.ir) then if(ra(j).lt.ra(j+1)) j=j + 1 endif if(rra.lt.ra(j)) then ra(i)=ra(j) i=j j=j+j else j=ir + 1 endif go to 20 endif ra(i)=rra go to 10 end c REFERENCE: Carnegie Mellon University Statistics Library. For c more information and access to subroutines, contact: c http://www.stat.cmu.edu/general/ c Additional references given in the prologue of each c subroutine. SUBROUTINE DQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB) C***BEGIN PROLOGUE DQRDC C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D5 C***KEYWORDS DECOMPOSITION,DOUBLE PRECISION,LINEAR ALGEBRA,LINPACK, C MATRIX,ORTHOGONAL TRIANGULAR C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE Uses Householder transformations to compute the Qr factori- C zation of N by P matrix X. Column pivoting is optional. C***DESCRIPTION C C DQRDC uses Householder transformations to compute the QR C factorization of an N by P matrix X. Column pivoting C based on the 2-norms of the reduced columns may be C performed at the user's option. C C On Entry C C X DOUBLE PRECISION(LDX,P), where LDX .GE. N. C X contains the matrix whose decomposition is to be C computed. C C LDX INTEGER. C LDX is the leading dimension of the array X. C C N INTEGER. C N is the number of rows of the matrix X. C C P INTEGER. C P is the number of columns of the matrix X. C C JPVT INTEGER(P). C JPVT contains integers that control the selection C of the pivot columns. The K-th column X(K) of X C is placed in one of three classes according to the C value of JPVT(K). C C If JPVT(K) .GT. 0, then X(K) is an initial C column. C C If JPVT(K) .EQ. 0, then X(K) is a free column. C C If JPVT(K) .LT. 0, then X(K) is a final column. C C Before the decomposition is computed, initial columns C are moved to the beginning of the array X and final C columns to the end. Both initial and final columns C are frozen in place during the computation and only C free columns are moved. At the K-th stage of the C reduction, if X(K) is occupied by a free column C it is interchanged with the free column of largest C reduced norm. JPVT is not referenced if C JOB .EQ. 0. C C WORK DOUBLE PRECISION(P). C WORK is a work array. WORK is not referenced if C JOB .EQ. 0. C C JOB INTEGER. C JOB is an integer that initiates column pivoting. C If JOB .EQ. 0, no pivoting is done. C If JOB .NE. 0, pivoting is done. C C On Return C C X X contains in its upper triangle the upper C triangular matrix R of the QR factorization. C Below its diagonal X contains information from C which the orthogonal part of the decomposition C can be recovered. Note that if pivoting has C been requested, the decomposition is not that C of the original matrix X but that of X C with its columns permuted as described by JPVT. C C QRAUX DOUBLE PRECISION(P). C QRAUX contains further information required to recover C the orthogonal part of the decomposition. C C JPVT JPVT(K) contains the index of the column of the C original matrix that has been interchanged into C the K-th column, if pivoting was requested. C C LINPACK. This version dated 08/14/78 . C G. W. Stewart, University of Maryland, Argonne National Lab. C C DQRDC uses the following functions and subprograms. C C BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2 C Fortran DABS,DMAX1,MIN0,DSQRT C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DDOT,DNRM2,DSCAL,DSWAP C***END PROLOGUE DQRDC INTEGER LDX,N,P,JOB INTEGER JPVT(1) DOUBLE PRECISION X(LDX,1),QRAUX(1),WORK(1) C INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU DOUBLE PRECISION MAXNRM,DNRM2,TT DOUBLE PRECISION DDOT,NRMXL,T LOGICAL NEGJ,SWAPJ C C***FIRST EXECUTABLE STATEMENT DQRDC PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 C C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. C DO 20 J = 1, P SWAPJ = JPVT(J) .GT. 0 NEGJ = JPVT(J) .LT. 0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J .NE. PL) CALL DSWAP(N,X(1,PL),1,X(1,J),1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ = 1, P J = P - JJ + 1 IF (JPVT(J) .GE. 0) GO TO 40 JPVT(J) = -JPVT(J) IF (J .EQ. PU) GO TO 30 CALL DSWAP(N,X(1,PU),1,X(1,J),1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE C C COMPUTE THE NORMS OF THE FREE COLUMNS. C IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUX(J) = DNRM2(N,X(1,J),1) WORK(J) = QRAUX(J) 70 CONTINUE 80 CONTINUE C C PERFORM THE HOUSEHOLDER REDUCTION OF X. C LUP = MIN0(N,P) DO 200 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 C C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. C MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (QRAUX(J) .LE. MAXNRM) GO TO 90 MAXNRM = QRAUX(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL DSWAP(N,X(1,L),1,X(1,MAXJ),1) QRAUX(MAXJ) = QRAUX(L) WORK(MAXJ) = WORK(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUX(L) = 0.0D0 IF (L .EQ. N) GO TO 190 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXL = DNRM2(N-L+1,X(L,L),1) IF (NRMXL .EQ. 0.0D0) GO TO 180 IF (X(L,L) .NE. 0.0D0) NRMXL = DSIGN(NRMXL,X(L,L)) CALL DSCAL(N-L+1,1.0D0/NRMXL,X(L,L),1) X(L,L) = 1.0D0 + X(L,L) C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. C LP1 = L + 1 IF (P .LT. LP1) GO TO 170 DO 160 J = LP1, P T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L) CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 150 IF (QRAUX(J) .EQ. 0.0D0) GO TO 150 TT = 1.0D0 - (DABS(X(L,J))/QRAUX(J))**2 TT = DMAX1(TT,0.0D0) T = TT TT = 1.0D0 + 0.05D0*TT*(QRAUX(J)/WORK(J))**2 IF (TT .EQ. 1.0D0) GO TO 130 QRAUX(J) = QRAUX(J)*DSQRT(T) GO TO 140 130 CONTINUE QRAUX(J) = DNRM2(N-L,X(L+1,J),1) WORK(J) = QRAUX(J) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE C C SAVE THE TRANSFORMATION. C QRAUX(L) = X(L,L) X(L,L) = -NRMXL 180 CONTINUE 190 CONTINUE 200 CONTINUE RETURN END SUBROUTINE DQRANK(A,LDA,M,N,TOL,KR,JPVT,QRAUX,WORK) C***BEGIN PROLOGUE DQRANK C***REVISION NOVEMBER 15, 1980 C***REFER TO DQRLSS C***KEYWORD(S) OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS C***AUTHOR D. KAHANER (NBS) C***DATE WRITTEN C***PURPOSE C COMPUTES THE QR FACTORIZATION OF AN M BY N MATRIX A. THIS C ROUTINE IS USED IN CONJUNCTION WITH DQRLSS TO SOLVE LINEAR C SYSTEMS OF EQUATIONS IN A LEAST SQUARE SENSE. C***DESCRIPTION C C DQRANK IS USED IN CONJUNCTION WITH DQRLSS TO SOLVE C OVERDETERMINED, UNDERDETERMINED AND SINGULAR LINEAR SYSTEMS C IN A LEAST SQUARES SENSE. C DQRANK USES THE LINPACK SUBROUTINE DQRDC TO COMPUTE THE QR C FACTORIZATION, WITH COLUMN PIVOTING, OF AN M BY N MATRIX A . C THE NUMERICAL RANK IS DETERMINED USING THE TOLERANCE TOL. C C FOR MORE INFORMATION, SEE CHAPTER 9 OF LINPACK USERS GUIDE, C J. DONGARRA ET ALL, SIAM, 1979. C C ON ENTRY C C A DOUBLE PRECISION (LDA,N) . C THE MATRIX WHOSE DECOMPOSITION IS TO BE COMPUTED. C C LDA INTEGER. C THE LEADING DIMENSION OF A . C C M INTEGER. C THE NUMBER OF ROWS OF A . C C N INTEGER. C THE NUMBER OF COLUMNS OF A . C C TOL DOUBLE PRECISION. C A RELATIVE TOLERANCE USED TO DETERMINE THE NUMERICAL C RANK. THE PROBLEM SHOULD BE SCALED SO THAT ALL THE ELEME C OF A HAVE ROUGHLY THE SAME ABSOLUTE ACCURACY, EPS. TH C REASONABLE VALUE FOR TOL IS ROUGHLY EPS DIVIDED BY C THE MAGNITUDE OF THE LARGEST ELEMENT. C C JPVT INTEGER(N) C C QRAUX DOUBLE PRECISION(N) C C WORK DOUBLE PRECISION(N) C THREE AUXILLIARY VECTORS USED BY DQRDC . C C ON RETURN C C A CONTAINS THE OUTPUT FROM DQRDC. C THE TRIANGULAR MATRIX R OF THE QR FACTORIZATION IS C CONTAINED IN THE UPPER TRIANGLE AND INFORMATION NEEDED C TO RECOVER THE ORTHOGONAL MATRIX Q IS STORED BELOW C THE DIAGONAL IN A AND IN THE VECTOR QRAUX . C C KR INTEGER. C THE NUMERICAL RANK. C C JPVT THE PIVOT INFORMATION FROM DQRDC. C C COLUMNS JPVT(1),...,JPVT(KR) OF THE ORIGINAL MATRIX ARE LINEARLY C INDEPENDENT TO WITHIN THE TOLERANCE TOL AND THE REMAINING COLUMNS C ARE LINEARLY DEPENDENT. ABS(A(1,1))/ABS(A(KR,KR)) IS AN ESTIMATE C OF THE CONDITION NUMBER OF THE MATRIX OF INDEPENDENT COLUMNS, C AND OF R . THIS ESTIMATE WILL BE .LE. 1/TOL . C C USAGE.....SEE SUBROUTINE DQRLSS C C***REFERENCE(S) C DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED DQRDC C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) DOUBLE PRECISION A(LDA,1),TOL,QRAUX(1),WORK(1) INTEGER J,K C***FIRST EXECUTABLE STATEMENT DO 10 J = 1, N JPVT(J) = 0 10 CONTINUE CALL DQRDC(A,LDA,M,N,QRAUX,JPVT,WORK,1) KR = 0 K = MIN0(M,N) DO 20 J = 1, K IF (ABS(A(J,J)) .LE. TOL*ABS(A(1,1))) GO TO 30 KR = J 20 CONTINUE 30 RETURN END SUBROUTINE DQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) C***BEGIN PROLOGUE DQRSL C***DATE WRITTEN 780814 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D9,D2A1 C***KEYWORDS DOUBLE PRECISION,LINEAR ALGEBRA,LINPACK,MATRIX, C ORTHOGONAL TRIANGULAR,SOLVE C***AUTHOR STEWART, G. W., (U. OF MARYLAND) C***PURPOSE Applies the output of DQRDC to compute coordinate C transformations, projections, and least squares solutions. C***DESCRIPTION C C DQRSL applies the output of DQRDC to compute coordinate C transformations, projections, and least squares solutions. C For K .LE. MIN(N,P), let XK be the matrix C C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C C formed from columnns JPVT(1), ... ,JPVT(K) of the original C N X P matrix X that was input to DQRDC (if no pivoting was C done, XK consists of the first K columns of X in their C original order). DQRDC produces a factored orthogonal matrix Q C and an upper triangular matrix R such that C C XK = Q * (R) C (0) C C This information is contained in coded form in the arrays C X and QRAUX. C C On Entry C C X DOUBLE PRECISION(LDX,P). C X contains the output of DQRDC. C C LDX INTEGER. C LDX is the leading dimension of the array X. C C N INTEGER. C N is the number of rows of the matrix XK. It must C have the same value as N in DQRDC. C C K INTEGER. C K is the number of columns of the matrix XK. K C must not be greater than MIN(N,P), where P is the C same as in the calling sequence to DQRDC. C C QRAUX DOUBLE PRECISION(P). C QRAUX contains the auxiliary output from DQRDC. C C Y DOUBLE PRECISION(N) C Y contains an N-vector that is to be manipulated C by DQRSL. C C JOB INTEGER. C JOB specifies what is to be computed. JOB has C the decimal expansion ABCDE, with the following C meaning. C C If A .NE. 0, compute QY. C If B,C,D, or E .NE. 0, compute QTY. C If C .NE. 0, compute B. C If D .NE. 0, compute RSD. C If E .NE. 0, compute XB. C C Note that a request to compute B, RSD, or XB C automatically triggers the computation of QTY, for C which an array must be provided in the calling C sequence. C C On Return C C QY DOUBLE PRECISION(N). C QY contains Q*Y, if its computation has been C requested. C C QTY DOUBLE PRECISION(N). C QTY contains TRANS(Q)*Y, if its computation has C been requested. Here TRANS(Q) is the C transpose of the matrix Q. C C B DOUBLE PRECISION(K) C B contains the solution of the least squares problem C C minimize norm2(Y - XK*B), C C if its computation has been requested. (Note that C if pivoting was requested in DQRDC, the J-th C component of B will be associated with column JPVT(J) C of the original matrix X that was input into DQRDC.) C C RSD DOUBLE PRECISION(N). C RSD contains the least squares residual Y - XK*B, C if its computation has been requested. RSD is C also the orthogonal projection of Y onto the C orthogonal complement of the column space of XK. C C XB DOUBLE PRECISION(N). C XB contains the least squares approximation XK*B, C if its computation has been requested. XB is also C the orthogonal projection of Y onto the column space C of X. C C INFO INTEGER. C INFO is zero unless the computation of B has C been requested and R is exactly singular. In C this case, INFO is the index of the first zero C diagonal element of R and B is left unaltered. C C The parameters QY, QTY, B, RSD, and XB are not referenced C if their computation is not requested and in this case C can be replaced by dummy variables in the calling program. C To save storage, the user may in some cases use the same C array for different parameters in the calling sequence. A C frequently occuring example is when one wishes to compute C any of B, RSD, or XB and does not need Y or QTY. In this C case one may identify Y, QTY, and one of B, RSD, or XB, while C providing separate arrays for anything else that is to be C computed. Thus the calling sequence C C CALL DQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) C C will result in the computation of B and RSD, with RSD C overwriting Y. More generally, each item in the following C list contains groups of permissible identifications for C a single calling sequence. C C 1. (Y,QTY,B) (RSD) (XB) (QY) C C 2. (Y,QTY,RSD) (B) (XB) (QY) C C 3. (Y,QTY,XB) (B) (RSD) (QY) C C 4. (Y,QY) (QTY,B) (RSD) (XB) C C 5. (Y,QY) (QTY,RSD) (B) (XB) C C 6. (Y,QY) (QTY,XB) (B) (RSD) C C In any group the value returned in the array allocated to C the group corresponds to the last member of the group. C C LINPACK. This version dated 08/14/78 . C G. W. Stewart, University of Maryland, Argonne National Lab. C C DQRSL uses the following functions and subprograms. C C BLAS DAXPY,DCOPY,DDOT C Fortran DABS,MIN0,MOD C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., C *LINPACK USERS GUIDE*, SIAM, 1979. C***ROUTINES CALLED DAXPY,DCOPY,DDOT C***END PROLOGUE DQRSL INTEGER LDX,N,K,JOB,INFO DOUBLE PRECISION X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1), 1 XB(1) C INTEGER I,J,JJ,JU,KP1 DOUBLE PRECISION DDOT,T,TEMP LOGICAL CB,CQY,CQTY,CR,CXB C C SET INFO FLAG. C C***FIRST EXECUTABLE STATEMENT DQRSL INFO = 0 C C DETERMINE WHAT IS TO BE COMPUTED. C CQY = JOB/10000 .NE. 0 CQTY = MOD(JOB,10000) .NE. 0 CB = MOD(JOB,1000)/100 .NE. 0 CR = MOD(JOB,100)/10 .NE. 0 CXB = MOD(JOB,10) .NE. 0 JU = MIN0(K,N-1) C C SPECIAL ACTION WHEN N=1. C IF (JU .NE. 0) GO TO 40 IF (CQY) QY(1) = Y(1) IF (CQTY) QTY(1) = Y(1) IF (CXB) XB(1) = Y(1) IF (.NOT.CB) GO TO 30 IF (X(1,1) .NE. 0.0D0) GO TO 10 INFO = 1 GO TO 20 10 CONTINUE B(1) = Y(1)/X(1,1) 20 CONTINUE 30 CONTINUE IF (CR) RSD(1) = 0.0D0 GO TO 250 40 CONTINUE C C SET UP TO COMPUTE QY OR QTY. C IF (CQY) CALL DCOPY(N,Y,1,QY,1) IF (CQTY) CALL DCOPY(N,Y,1,QTY,1) IF (.NOT.CQY) GO TO 70 C C COMPUTE QY. C DO 60 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 50 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QY(J),1) X(J,J) = TEMP 50 CONTINUE 60 CONTINUE 70 CONTINUE IF (.NOT.CQTY) GO TO 100 C C COMPUTE TRANS(Q)*Y. C DO 90 J = 1, JU IF (QRAUX(J) .EQ. 0.0D0) GO TO 80 TEMP = X(J,J) X(J,J) = QRAUX(J) T = -DDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,QTY(J),1) X(J,J) = TEMP 80 CONTINUE 90 CONTINUE 100 CONTINUE C C SET UP TO COMPUTE B, RSD, OR XB. C IF (CB) CALL DCOPY(K,QTY,1,B,1) KP1 = K + 1 IF (CXB) CALL DCOPY(K,QTY,1,XB,1) IF (CR .AND. K .LT. N) CALL DCOPY(N-K,QTY(KP1),1,RSD(KP1),1) IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120 DO 110 I = KP1, N XB(I) = 0.0D0 110 CONTINUE 120 CONTINUE IF (.NOT.CR) GO TO 140 DO 130 I = 1, K RSD(I) = 0.0D0 130 CONTINUE 140 CONTINUE IF (.NOT.CB) GO TO 190 C C COMPUTE B. C DO 170 JJ = 1, K J = K - JJ + 1 IF (X(J,J) .NE. 0.0D0) GO TO 150 INFO = J C ......EXIT GO TO 180 150 CONTINUE B(J) = B(J)/X(J,J) IF (J .EQ. 1) GO TO 160 T = -B(J) CALL DAXPY(J-1,T,X(1,J),1,B,1) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 C C COMPUTE RSD OR XB AS REQUIRED. C DO 230 JJ = 1, JU J = JU - JJ + 1 IF (QRAUX(J) .EQ. 0.0D0) GO TO 220 TEMP = X(J,J) X(J,J) = QRAUX(J) IF (.NOT.CR) GO TO 200 T = -DDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 200 CONTINUE IF (.NOT.CXB) GO TO 210 T = -DDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) CALL DAXPY(N-J+1,T,X(J,J),1,XB(J),1) 210 CONTINUE X(J,J) = TEMP 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE RETURN END SUBROUTINE DQRLSS(A,LDA,M,N,KR,B,X,RSD,JPVT,QRAUX) C***BEGIN PROLOGUE DQRLSS C***REVISION NOVEMBER 15, 1980 C***CATEGORY NO. D9 C***KEYWORD(S) OVERDETERMINED,LEAST SQUARES,LINEAR EQUATIONS C***AUTHOR D. KAHANER (NBS) C***DATE WRITTEN C***PURPOSE C SOLVES AN OVERDETERMINED, UNDERDETERMINED, OR SINGULAR SYSTEM OF C LINEAR EQUATIONS IN LEAST SQUARE SENSE. C***DESCRIPTION C C DQRLSS USES THE LINPACK SUBROUTINE DQRSL TO SOLVE AN OVERDETERMINE C UNDERDETERMINED, OR SINGULAR LINEAR SYSTEM IN A LEAST SQUARES C SENSE. IT MUST BE PRECEEDED BY DQRANK . C THE SYSTEM IS A*X APPROXIMATES B WHERE A IS M BY N WITH C RANK KR (DETERMINED BY DQRANK), B IS A GIVEN M-VECTOR, C AND X IS THE N-VECTOR TO BE COMPUTED. A SOLUTION X WITH C AT MOST KR NONZERO COMPONENTS IS FOUND WHICH MINIMIMZES C THE 2-NORM OF THE RESIDUAL, A*X - B . C C ON ENTRY C C A,LDA,M,N,KR,JPVT,QRAUX C THE OUTPUT FROM DQRANK . C C B DOUBLE PRECISION(M) . C THE RIGHT HAND SIDE OF THE LINEAR SYSTEM. C C ON RETURN C C X DOUBLE PRECISION(N) . C A LEAST SQUARES SOLUTION TO THE LINEAR SYSTEM. C C RSD DOUBLE PRECISION(M) . C THE RESIDUAL, B - A*X . RSD MAY OVERWITE B . C C USAGE.... C ONCE THE MATRIX A HAS BEEN FORMED, DQRANK SHOULD BE C CALLED ONCE TO DECOMPOSE IT. THEN FOR EACH RIGHT HAND C SIDE, B, DQRLSS SHOULD BE CALLED ONCE TO OBTAIN THE C SOLUTION AND RESIDUAL. C C C***REFERENCE(S) C DONGARRA, ET AL, LINPACK USERS GUIDE, SIAM, 1979 C***ROUTINES CALLED DQRSL C***END PROLOGUE INTEGER LDA,M,N,KR,JPVT(1) INTEGER J,K,INFO DOUBLE PRECISION T DOUBLE PRECISION A(LDA,1),B(1),X(1),RSD(1),QRAUX(1) C***FIRST EXECUTABLE STATEMENT IF (KR .NE. 0) 1 CALL DQRSL(A,LDA,M,KR,QRAUX,B,RSD,RSD,X,RSD,RSD,110,INFO) DO 40 J = 1, N JPVT(J) = -JPVT(J) IF (J .GT. KR) X(J) = 0.0 40 CONTINUE DO 70 J = 1, N IF (JPVT(J) .GT. 0) GO TO 70 K = -JPVT(J) JPVT(J) = K 50 CONTINUE IF (K .EQ. J) GO TO 60 T = X(J) X(J) = X(K) X(K) = T JPVT(K) = -JPVT(K) K = JPVT(K) GO TO 50 60 CONTINUE 70 CONTINUE RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DAXPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A7 C***KEYWORDS BLAS,DOUBLE PRECISION,LINEAR ALGEBRA,TRIAD,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P computation y = a*x + y C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C DA double precision scalar multiplier C DX double precision vector with N elements C INCX storage spacing between elements of DX C DY double precision vector with N elements C INCY storage spacing between elements of DY C C --Output-- C DY double precision result (unchanged if N .LE. 0) C C Overwrite double precision DY with double precision DA*DX + DY. C For I = 0 to N-1, replace DY(LY+I*INCY) with DA*DX(LX+I*INCX) + C DY(LY+I*INCY), where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N C and LY is defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DAXPY C DOUBLE PRECISION DX(1),DY(1),DA C***FIRST EXECUTABLE STATEMENT DAXPY IF(N.LE.0.OR.DA.EQ.0.D0) RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. 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 C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. 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 C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX DY(I) = DA*DX(I) + DY(I) 70 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DCOPY C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A5 C***KEYWORDS BLAS,COPY,DOUBLE PRECISION,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. vector copy y = x C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C DX double precision vector with N elements C INCX storage spacing between elements of DX C DY double precision vector with N elements C INCY storage spacing between elements of DY C C --Output-- C DY copy of vector DX (unchanged if N .LE. 0) C C Copy double precision DX to double precision DY. C For I = 0 to N-1, copy DX(LX+I*INCX) to DY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DCOPY C DOUBLE PRECISION DX(1),DY(1) C***FIRST EXECUTABLE STATEMENT DCOPY IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 7. C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 60 CONTINUE NS=N*INCX DO 70 I=1,NS,INCX DY(I) = DX(I) 70 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C***BEGIN PROLOGUE DSCAL C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A6 C***KEYWORDS BLAS,LINEAR ALGEBRA,SCALE,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. vector scale x = a*x C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C DA double precision scale factor C DX double precision vector with N elements C INCX storage spacing between elements of DX C C --Output-- C DX double precision result (unchanged if N.LE.0) C C Replace double precision DX by double precision DA*DX. C For I = 0 to N-1, replace DX(1+I*INCX) with DA * DX(1+I*INCX) C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DSCAL C DOUBLE PRECISION DA,DX(1) C***FIRST EXECUTABLE STATEMENT DSCAL IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENTS NOT EQUAL TO 1. C NS = N*INCX DO 10 I = 1,NS,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. 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 SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DSWAP C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A5 C***KEYWORDS BLAS,DOUBLE PRECISION,INTERCHANGE,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE Interchange d.p. vectors C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C DX double precision vector with N elements C INCX storage spacing between elements of DX C DY double precision vector with N elements C INCY storage spacing between elements of DY C C --Output-- C DX input vector DY (unchanged if N .LE. 0) C DY input vector DX (unchanged if N .LE. 0) C C Interchange double precision DX and double precision DY. C For I = 0 to N-1, interchange DX(LX+I*INCX) and DY(LY+I*INCY), C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DSWAP C DOUBLE PRECISION DX(1),DY(1),DTEMP1,DTEMP2,DTEMP3 C***FIRST EXECUTABLE STATEMENT DSWAP IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 DTEMP1 = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP1 IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3. C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP1 = DX(I) DTEMP2 = DX(I+1) DTEMP3 = DX(I+2) DX(I) = DY(I) DX(I+1) = DY(I+1) DX(I+2) = DY(I+2) DY(I) = DTEMP1 DY(I+1) = DTEMP2 DY(I+2) = DTEMP3 50 CONTINUE RETURN 60 CONTINUE C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C NS = N*INCX DO 70 I=1,NS,INCX DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C***BEGIN PROLOGUE DDOT C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A4 C***KEYWORDS BLAS,DOUBLE PRECISION,INNER PRODUCT,LINEAR ALGEBRA,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE D.P. inner product of d.p. vectors C***DESCRIPTION C C B L A S Subprogram C Description of Parameters C C --Input-- C N number of elements in input vector(s) C DX double precision vector with N elements C INCX storage spacing between elements of DX C DY double precision vector with N elements C INCY storage spacing between elements of DY C C --Output-- C DDOT double precision dot product (zero if N .LE. 0) C C Returns the dot product of double precision DX and DY. C DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY) C where LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is C defined in a similar way using INCY. C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DDOT C DOUBLE PRECISION DX(1),DY(1) C***FIRST EXECUTABLE STATEMENT DDOT DDOT = 0.D0 IF(N.LE.0)RETURN IF(INCX.EQ.INCY) IF(INCX-1) 5,20,60 5 CONTINUE C C CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS. 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 DDOT = DDOT + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1. C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DDOT = DDOT + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + 1 DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 60 CONTINUE NS = N*INCX DO 70 I=1,NS,INCX DDOT = DDOT + DX(I)*DY(I) 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DNRM2(N,DX,INCX) C***BEGIN PROLOGUE DNRM2 C***DATE WRITTEN 791001 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. D1A3B C***KEYWORDS BLAS,DOUBLE PRECISION,EUCLIDEAN,L2,LENGTH,LINEAR ALGEBRA, C NORM,VECTOR C***AUTHOR LAWSON, C. L., (JPL) C HANSON, R. J., (SNLA) C KINCAID, D. R., (U. OF TEXAS) C KROGH, F. T., (JPL) C***PURPOSE Euclidean length (L2 norm) of d.p. vector C***DESCRIPTION C C B L A S Subprogram C Description of parameters C C --Input-- C N number of elements in input vector(s) C DX double precision vector with N elements C INCX storage spacing between elements of DX C C --Output-- C DNRM2 double precision result (zero if N .LE. 0) C C Euclidean norm of the N-vector stored in DX() with storage C increment INCX . C If N .LE. 0 return with result = 0. C If N .GE. 1 then INCX must be .GE. 1 C C C.L. Lawson, 1978 Jan 08 C C Four phase method using two built-in constants that are C hopefully applicable to all machines. C CUTLO = maximum of DSQRT(U/EPS) over all known machines. C CUTHI = minimum of DSQRT(V) over all known machines. C where C EPS = smallest no. such that EPS + 1. .GT. 1. C U = smallest positive no. (underflow limit) C V = largest no. (overflow limit) C C Brief outline of algorithm.. C C Phase 1 scans zero components. C move to phase 2 when a component is nonzero and .LE. CUTLO C move to phase 3 when a component is .GT. CUTLO C move to phase 4 when a component is .GE. CUTHI/M C where M = N for X() real and M = 2*N for complex. C C Values for CUTLO and CUTHI.. C From the environmental parameters listed in the IMSL converter C document the limiting values are as followS.. C CUTLO, S.P. U/EPS = 2**(-102) for Honeywell. Close seconds are C Univac and DEC at 2**(-103) C Thus CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 for Univac, Honeywell, and DEC. C Thus CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) for Honeywell and DEC. C Thus CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. same as S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C***REFERENCES LAWSON C.L., HANSON R.J., KINCAID D.R., KROGH F.T., C *BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE*, C ALGORITHM NO. 539, TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOLUME 5, NUMBER 3, SEPTEMBER 1979, 308-323 C***ROUTINES CALLED (NONE) C***END PROLOGUE DNRM2 INTEGER NEXT DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C***FIRST EXECUTABLE STATEMENT DNRM2 IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END REAL FUNCTION UNI(JD) C***BEGIN PROLOGUE UNI C***DATE WRITTEN 810915 C***REVISION DATE 830805 C***CATEGORY NO. L6A21 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS C***AUTHOR BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS C KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV C C***PURPOSE THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1 C AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS C AT LEAST AS LARGE AS 32767. C***DESCRIPTION C C THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER C [0,1). IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS C INTEGERS AT LEAST AS LARGE AS 32767. C C C USE C FIRST TIME.... C Z = UNI(JD) C HERE JD IS ANY N O N - Z E R O INTEGER. C THIS CAUSES INITIALIZATION OF THE PROGRAM C AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z. C SUBSEQUENT TIMES... C Z = UNI(0) C CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z. C C C.................................................................. C NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER C TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION..... C C MACHINE DEPENDENCIES... C MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE C FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT. C THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED C IN LINE WITH REMARK A BELOW. C C REMARKS... C A. THIS PROGRAM CAN BE USED IN TWO WAYS: C (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS, C SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR, C (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE C GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE C LARGEST POSSIBLE VALUE. C B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL C INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'. C IF MDIG=16 ONE SHOULD FIND THAT C THE FIRST EVALUATION C Z=UNI(305) GIVES Z=.027832881... C THE SECOND EVALUATION C Z=UNI(0) GIVES Z=.56102176... C THE THIRD EVALUATION C Z=UNI(0) GIVES Z=.41456343... C THE THOUSANDTH EVALUATION C Z=UNI(0) GIVES Z=.19797357... C C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. C***ROUTINES CALLED I1MACH,XERROR C***END PROLOGUE UNI INTEGER M(17) C SAVE I,J,M,M1,M2 C DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11), 1 M(12),M(13),M(14),M(15),M(16),M(17) 2 / 30788,23052,2053,19346,10646,19427,23975, 3 19049,10949,19693,29746,26748,2796,23890, 4 29168,31924,16499 / DATA M1,M2,I,J / 32767,256,5,17 / C***FIRST EXECUTABLE STATEMENT UNI IF(JD .EQ. 0) GO TO 3 C FILL MDIG=16 C BE SURE THAT MDIG AT LEAST 16... M1= 2**(MDIG-2) + (2**(MDIG-2)-1) M2 = 2**(MDIG/2) JSEED = MIN0(IABS(JD),M1) IF( MOD(JSEED,2).EQ.0 ) JSEED=JSEED-1 K0 =MOD(9069,M2) K1 = 9069/M2 J0 = MOD(JSEED,M2) J1 = JSEED/M2 DO 2 I=1,17 JSEED = J0*K0 J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2) J0 = MOD(JSEED,M2) 2 M(I) = J0+M2*J1 I=5 J=17 C BEGIN MAIN LOOP HERE 3 K=M(I)-M(J) IF(K .LT. 0) K=K+M1 M(J)=K I=I-1 IF(I .EQ. 0) I=17 J=J-1 IF(J .EQ. 0) J=17 UNI=FLOAT(K)/FLOAT(M1) RETURN END c ******************************************************** Character*64 Function TRANLC(text) C C PURPOSE: Tranlate a character variable to all lowercase letters. C Non-alphabetic characters are not affected. C C COMMENTS: Text = lowercase version of itself C TRANLC = original version of text C C CODE DEPENDENCIES: C Routine Name File C N/A C C AUTHOR: Robert D. Stewart C DATE: May 19, 1992 C Modified by: JJP 1-21-99 C CHARACTER*(*) text Character*64 text2 INTEGER iasc,i,ilen text2 = text ilen = LEN(text) DO I=1,ilen iasc = ICHAR(text(i:i)) IF ((iasc.GT.64).AND.(iasc.LT.91)) THEN text(i:i) = CHAR(iasc+32) ELSE text(i:i) = CHAR(iasc) ENDIF ENDDO TRANLC = text2(1:ilen) RETURN END c ******************************************************** c ******************************************************** SUBROUTINE SQUEEZE(text) C C PURPOSE: Remove any extra and leading blanks from character C variable text. C C COMMENTS: The original "text" is destroyed on exit. This C routine is superseded by the RmvRep routines, but C is included for backwards compatibility. C C CODE DEPENDENCIES: C Routine Name File C N/A C C AUTHOR: Robert D. Stewart C DATE: May 19, 1992 C Modified: JJP 1-21-98 because it did not work. removed leading C blanks erasing thing C CHARACTER*(*) text Character*64 text2 CHARACTER*1 ch INTEGER iasc,icnt,ilen,i if(text(1:1).eq.' ')RETURN ilen = LEN(text) IF (ilen.GT.1) THEN icnt = 1 DO I=1,ilen iasc = ICHAR(text(i:i)) if (iasc.ne.32) then ch = text(i:i) text2(icnt:icnt) = ch icnt = icnt + 1 ENDIF ENDDO DO I=icnt,ilen text2(i:i) = ' ' ENDDO text = '' ilen = LEN(text2) text(1:ilen) = text2(1:ilen) ENDIF RETURN END