Richard, The software (HYCHOL.F) was written by H. Klein almost 20 years ago. I got the code from N. Zacharias (USNO), who apparently told you that I had tested a version of HYCHOL. I have a brief manual for HYCHOL, written in German. For our purposes, I translated it into English; I am including the file (hychol.tex) and will send you a FAX of the original, which also contains some useful diagrams. In addition, I am including a small test program. It works on a hard- wired matrix, which is included in the FAX. Files: HYCHOL.F: Subroutines MATOUT.F: additional subroutine, apparently written by N. Zacharias; unimportant MAIN.F: N. Zacharias' main program for plate reduction; sets up matrices TEST.F small test program: should print out: 1 2 3 4 5 6 7 8 (approximately) All the programs ran on my 386 PC, under NDP Fortran 1.4e (ancient). I never used it for serious data reduction. I am not using this software anymore, because I have developed my own C-code, which does a Cholesky decomposition of a sparse matrix given in row-wise format. This program, however, is in an experimental stage. The algorithm is quite simple, though. Greetings, Christian Hummel US Naval Observatory - AD5 3450 Massachusetts Avenue NW Washington DC 20392-5420 Tel.: (202) 653 0950 FAX : (202) 653 0941 e-mail: cah@fornax.usno.navy.mil ---------------------------------------HYCHOL.TEX (LaTex) \documentstyle[12pt]{article} \textwidth 15cm \textheight 24cm \topmargin -2cm \oddsidemargin -0cm \evensidemargin -0cm \begin{document} \begin{center} {\Large\bf User's guide to HYCHOL} \end{center} \section{Purpose of the program} The program solves large, symmetric, and positive definite sets of linear equations with non-singular design matrix. The design matrix does not have to fit completely into the computer memory. Coefficients are stored on disk with special I/O routines. More than one right-hand side can be specified. By using unit vectors as right-hand side, the program can be used to partially invert the matrix. The program is optimized for matrices with band structure. The source code is written in FORTRAN. \section{Mathematical procedure} The design matrix of a set of linear equations can be written as a hyper-matrix of sub-matrices, using a symmetric grid. The system of sub-matrices is dealt with a modified Cholesky-algorithm. The original matrix, \begin{equation} \sum_{k=1}^n A_{ik} X_k = A_i,\;\;\;i=1,n, \end{equation} is transformed into \begin{equation} \sum_{k=i}^n R_{ik} X_k = R_i,\;\;\;i=1,n,\;k=i,n, \end{equation} with \[\left. \begin{array}{r} A_{ik}=\\A_i=\\R_{ik}=\\R_i=\end{array}\right\} {\rm sub-matrices\;of\;the} \left\{ \begin{array}{c} {\rm left}\\{\rm right}\\{\rm left}\\{\rm right} \end{array} \right\} {\rm sides\;of\;the} \left\{ \begin{array}{c}{\rm original}\\{\rm re-ordered}\end{array} \right\} {\rm system} \] \[X_k={\rm sub-matrices\;of\;the\;solution\;matrix\;of\;the\;set \;of\;linear\;equations}.\] $A_i$, $R_i$, and $X_k$ reduce to simple vectors in case there is only one right-hand side. The transformation of Eq.\ (2) into (1) is done with the following algorithm (omitting right-hand side index $k$): \begin{equation} R_{ik}=R_{ii}^{{\rm T}^{-1}}\cdot\left(A_{ik}-\sum_{r=1}^{i-1} R_{ri}^{\rm T}\cdot R_{rk} \right),\;\;\;i=1,n\;k=i,n. \end{equation} According to the aforementioned conditions, the design matrix is positive definite. Thus, all determinants of the main sections are positive. Since these do not change through the transformations of the elimination process, all diagonal sub-matrices are also positive definite. The Cholesky-decomposition of the latter can therefore be performed. The unknown sub-matrices $X_k$ can be calculated recursively from Eq.\ (2), starting with $X_n$. \section{Organization} The size of the sub-matrices can be chosen arbitrarily. Not more than three sub-matrices will be stored in memory at the same time. Thus, even very large sets of linear equations can be solved with limited computer memory (allowing for more I/O time, of course). The diagonal form of the sub-matrices $R_{ii}$ and the symmetry of the matrix-products in Eq.\ (3) for $k=i$ is accounted for in all matrix operations. Unlike indicated in Eqs.\ (1-3), the creation of indices for the sub-matrices follows common schemes for symmetric band matrices, assigning column index 1 to the diagonal matrices and continuing the count to the right. Corresponding to the structure of the hyper-matrix, an address-matrix has to be created, so that each sub-matrix is assigned a place in the address-matrix. An I/O-subroutine package handles all organizational matters of the matrix solution and the data transfer from and to the disk. \section{How to call HYCHOL and the meaning of the parameters} \noindent CALL HYCHOL (IADM, IHZEI, IHBA1, MAXZEI, IHBDI, IHBRA, A, B, C, IDIM,\\ \hspace*{72pt} IR, IGL, IGR) \begin{tabular}{ll} IADM (IHZEI, IHBA1) & address-matrix.\\ MAXZEI, IHBDI, IHBRA & maximum index value in left side \\ & of address-matrix IADM (I, K),\\ & I$=$1,..,MAXZEI;\\& K$=$1,...,IHBA1$-$1.\\ & For right side, I$=$1,...,MAXZEI; \\ & K$=$IHBA1.\\ A(IDIM,IDIM), B(,), C(,) & sub-matrices.\\ IR & maximum number of right-hand sides.\\ IGL, IGR & logical units for files of left and right sides. \end{tabular} Conditions: MAXZEI$\leq$IHZEI; IHBRA$+$IHBDI$\leq$IHBA1$-$1; IR$\leq$IDIM. Before HYCHOL is called, the sub-matrices of the left and right-hand sides have to be stored to files IGL and IGR using subroutine PUTAR. It is not necessary to store null matrices. \section{The I/O-package} In order to handle the data transfer from and to the external storage and to organize the solution, an I/O-subroutine package is used. For this, the following is needed: \begin{itemize} \item three labeled common blockes:\\ INTEGER NOFIL,IFILE(7,2),IOU,KKKK,LPRU\\ COMMON/DEFIE/NOFIL,IFILE\\ COMMON/OUT/IOU,KKKK\\ COMMON/IBMF/LPRU \item assignments and other statements:\\ CHARACTER*8 LFNAM(2)\\ INTEGER LFILE(3,2)\\ IOU=6\\ IGL=17\\ IGR=18\\ OPEN(UNIT=IGL,ACCESS='DIRECT',RECL=19040)\\ OPEN(UNIT=IGR,ACCESS='DIRECT',RECL=19040)\\ CALL DINIT(2)\\ LFNAM(1)='LE.SCR'\\ LFILE(1,1)=IGL\\ LFILE(2,1)=IDIM*IDIM\\ LFILE(3,1)=IHZEI*(IHBA1-1)\\ LFNAM(2)='RI.SCR'\\ LFILE(1,2)=IGR\\ LFILE(2,2)=IDIM\\ LFILE(3,2)=IHZEI\\ CALL VORHY(IADM,IHZEI,IHBA1,LFILE,LFNAM) \item storing the sub-matrices of the left and right side:\\ CALL PUTAR(IADM,IHZEI,I,K,A,IDIM,IR,IA,IG)\\ \begin{tabular}{ll} I, K: & Address of submatrix in hyper-matrix IADM \\ IDIM, IR: & Dimension of submatrix A, IR$=$IDIM for\\ & sub-matrices of the left side\\ IA: & valid number of columns in A (IA$\leq$IR)\\ IG: & Logical file number (IGL/IGR) \end{tabular} \item reading of the results after solution with HYCHOL:\\ CALL GETAR(IADM,IHZEI,I,IHBA1,A,IDIM,IR,IGR) \item noting that:\\ The files for the left and right side of the set of equations are random direct access files, which are initialized by DINIT(2). They are opened by DOPEN(IUNIT,3HLFN,LRECL,MAXRE,MAXWRI) as sequential files, but are immediately transformed into random files, with\\ \begin{tabular}{ll} IUNIT$=$ & logical file number\\ LFN$=$ & logical file name\\ LRECL$=$ & record length in 60-bit words or double precision\\ MAXRE$=$ & maximum number of records\\ MAXWRI$=$ & number of words already written to file.\\ \end{tabular} MAXWRI is zero initially. DOPEN is executed in HYCHOL by calling VORHY, passing array LFILE. DEVICT(IUNIT) can be used to close a file, keeping the logical file name. Reading and writing to the direct access files is done using GETR(IUNIT,NRECD,FWA) and PUTR(IUNIT,NRECD,FWA), where NRECD$=$record number, FWA$=$first word address of the field being written or read. Records must have constant length. \end{itemize} \end{document} ----------------------------------TEST.F program test c c Program to test the HYCHOL subroutines with hard-wired matrix c c Declarations for PUTAR c c Address matrix: integer*4 ihzei,ihba1 parameter (ihzei=4,ihba1=5) integer*4 iadm(ihzei,ihba1) c Addresses: integer*4 iad,jad c Submatrix: integer*4 idim,ir parameter (idim=2,ir=1) real*8 ar(idim,ir) c Logical files and common blocks: integer*4 iou,igl,igr integer*4 kkkk,nofil,lpru,ifile(7,2),lfile(3,2) common/out/iou,kkkk common/ibmf/lpru common/defie/nofil,ifile character*8 lfnam(2) c c Declarations for HCHOL c c Size of hyper matrix, # sub.matr. in row, # band cols.: integer*4 maxzei,ihbdi,ihbra parameter (maxzei=4,ihbdi=1,ihbra=1) c Number of right hand sides: integer*4 iresei parameter (iresei=1) c Submatrices (work space): real*8 a(idim,idim),b(idim,idim),c(idim,idim) c c Function ixy10: integer*4 ixy10 external ixy10 c iou=6 igl=17 igr=18 ntrsz=32258 lpru=ntrsz/8 OPEN (igl,ACCESS='DIRECT',RECL=19040) OPEN (igr,ACCESS='DIRECT',RECL=19040) c c Initialize file tables call dinit(2) lfnam(1)='LE.SCR' lfile(1,1)=igl lfile(2,1)=idim*idim lfile(3,1)=ihzei*(ihba1-1) lfnam(2)='RI.SCR' lfile(1,2)=igr lfile(2,2)=idim lfile(3,2)=ihzei call vorhy(iadm,ihzei,ihba1,lfile,lfnam) c c Store sumatrices c submatrix 1,1 a(1,1)=0.24d0 a(2,1)=0.44d0 a(1,2)=0.44d0 a(2,2)=1.15d0 ia=idim call putar(iadm,ihzei,1,1,a,idim,idim,ia,igl) c Submatrix 2,1 a(1,1)=0.75d0 a(2,1)=0.83d0 a(1,2)=0.83d0 a(2,2)=1.58d0 call putar(iadm,ihzei,2,1,a,idim,idim,ia,igl) c Submatrix 3,1 a(1,1)=1.36d0 a(2,1)=0.73d0 a(1,2)=0.73d0 a(2,2)=1.40d0 call putar(iadm,ihzei,3,1,a,idim,idim,ia,igl) c border band a(1,1)=0.16d0 a(2,1)=0.47d0 a(1,2)=0.32d0 a(2,2)=0.61d0 call putar(iadm,ihzei,1,ihba1-1,a,idim,idim,ia,igl) a(1,1)=0.27d0 a(2,1)=0.57d0 a(1,2)=0.56d0 a(2,2)=0.72d0 call putar(iadm,ihzei,2,ihba1-1,a,idim,idim,ia,igl) a(1,1)=1.19d0 a(2,1)=1.66d0 a(1,2)=1.45d0 a(2,2)=1.18d0 call putar(iadm,ihzei,3,ihba1-1,a,idim,idim,ia,igl) c lower right a(1,1)=2.68d0 a(2,1)=2.33d0 a(1,2)=2.33d0 a(2,2)=3.08d0 call putar(iadm,ihzei,ihzei,ihba1-1,a,idim,idim,ia,igl) c c Store right hand side ia=ir ar(1,1)= 4.80d0 ar(2,1)=10.91d0 call putar(iadm,ihzei,1,ihba1,ar,idim,ir,ia,igr) ar(1,1)=11.94d0 ar(2,1)=18.56d0 call putar(iadm,ihzei,2,ihba1,ar,idim,ir,ia,igr) ar(1,1)=31.11d0 ar(2,1)=33.11d0 call putar(iadm,ihzei,3,ihba1,ar,idim,ir,ia,igr) c lower right ar(1,1)=57.50d0 ar(2,1)=61.38d0 call putar(iadm,ihzei,ihzei,ihba1,ar,idim,ir,ia,igr) c c call HYCHOL routines call hchol(iadm,ihzei,ihba1,maxzei,ihbdi,ihbra,a,b,c,idim, & iresei,igl,igr) c c get results do i=1,maxzei call getar(iadm,ihzei,i,ihba1,ar,idim,ir,igr) do j=1,idim write(6,*)ar(j,ir) enddo enddo c end ------------------------------------------------HYCHOL.F C H.B.A.P.P. --S80HYC-- MAIN PGM = LAST UPDATE: 21.10.87 C SUBROUTINE PACKAGE TO SOLVE LARGE SYSTEMS OF LINEAR EQUATIONS C WITH HYPER CHOLESKY METHODS. C SUBROUTINE GETR (IG,IRNU,A,I1,I2) C C READ RECORD IRNU OF FILE IG INTO A FOR DIO PACKAGE C DOUBLE PRECISION A,TR COMMON/DEFIE/NOFIL,IFILE(7,2) COMMON/OUT/IOU,KKKK DIMENSION A(I1),TR(2380) DIMENSION IA(3) CALL SEARC(IG,IPOS) IF(IPOS.EQ.0) CALL DERRO(1,IG,IA) K=I1*I2 IA(1)=IRNU IF(IRNU.GT.IFILE(7,IPOS)) CALL DERRO(4,IG,IA) IF(K.EQ.IFILE(3,IPOS)) GOTO 10 IA(2)=K IA(3)=IFILE(3,IPOS) CALL DERRO(5,IG,IA) 10 CONTINUE K1=1+(IRNU-1)*IFILE(5,IPOS) K2=K1+IFILE(5,IPOS)-1 L1=1 DO 2 M=K1,K2 L2=L1+IFILE(6,IPOS)-1 IF(L2.GT.K) L2=K READ(IG,REC=M) TR IND=0 DO 3 L=L1,L2 IND=IND+1 3 A(L)=TR(IND) 2 L1=L2+1 RETURN END SUBROUTINE PUTR(IG,IRNU,A,I1,I2) C C WRITE RECORD IRNU FROM A TO FILE IG FOR DIO-PACKAGE C DOUBLE PRECISION A,TR DIMENSION A(I1),TR(2380) COMMON/DEFIE/NOFIL,IFILE(7,2) COMMON/OUT/IOU,KKKK DIMENSION IA(3) CALL SEARC(IG,IPOS) IF(IPOS.EQ.0) CALL DERRO(1,IG,IA) K=I1*I2 IA(1)=IRNU IF(IRNU.GT.IFILE(4,IPOS)) CALL DERRO(7,IG,IA) IF(IRNU.GT.IFILE(7,IPOS)) IFILE(7,IPOS)=IRNU IF(K.EQ.IFILE(3,IPOS)) GOTO 10 IA(2)=K IA(3)=IFILE(3,IPOS) CALL DERRO(6,IG,IA) 10 CONTINUE K1=1+(IRNU-1)*IFILE(5,IPOS) K2=K1+IFILE(5,IPOS)-1 L1=1 DO 2 M=K1,K2 L2=L1+IFILE(6,IPOS)-1 IF(L2.GT.K) L2=K IND=0 DO 3 L=L1,L2 IND=IND+1 3 TR(IND)=A(L) WRITE(IG,REC=M) TR 2 L1=L2+1 RETURN END SUBROUTINE DOPEN(IUNIT,LFN,LRECL,MAXREC,MWREC) C C CREATE FILE TABLE FOR FILE IUNIT FOR DIO-PACKAGE C CHARACTER*8 LFN COMMON/DEFIE/NOFIL,IFILE(7,2) COMMON/OUT/IOU,KKKK COMMON/IBMF/LPRU CALL SEARC(IUNIT,IPOS) IF(IPOS.NE.0) CALL DERRO(3,IUNIT,NR) CALL SEARC(0,IPOS) IF(IPOS.EQ.0) CALL DERRO(2,IUNIT,NR) IFILE(1,IPOS)=IUNIT C IFILE(2,IPOS)=LFN IFILE(3,IPOS)=LRECL IFILE(5,IPOS)=(LRECL-1)/LPRU+1 IFILE(7,IPOS)=MWREC IFILE(4,IPOS)=MAXREC IFILE(6,IPOS)=LPRU C OPEN(UNIT=IUNIT,ERR=100,STATUS='SCRATCH',ACCESS='DIRECT', C 1IOSTAT=IOS,RECL=8*LRECL) RETURN 100 WRITE(IOU,200) IUNIT,LFN,IOS 200 FORMAT(/,27H >>>>> OPEN ERROR FOR FILE ,I4,3H, ,A8,3H, , 110H IOSTAT= ,I5,6H <<<<<) STOP END C SUBROUTINE HCHOL(IADM,N,MM1,NN,MB,MR,A,B,C,ID,IR,IGL,IGR) C C HYPERCHOLESKI METHOD FOR SOLUTION OF POSITIVE DEFINIT LINEAR C EQUATION SYSTEM C DOUBLE PRECISION A,B,C COMMON/OUT/IOU,KKKK DIMENSION IADM(N,MM1) DIMENSION A(ID,ID),B(ID,ID),C(ID,ID) KKKK=0 DO 21 KZ=1,NN K=KZ IF(KZ.GT.NN-MR) K=N-NN+KZ KS=1 IF(KZ.GT.NN-MR) KS=MM1-1+KZ-NN I1=IXY10(IADM,N,K,KS) CALL GETAR(IADM,N,K,KS,A,ID,ID,IGL) CALL XY9(A,I1,ID) IF(KKKK.NE.0) THEN PRINT*,'-- IN , AFTER , K,KS = ',K,KS PRINT*,'-- (CHOLESKY PARTITIONING) A(I,I) =' WRITE (6,'(1X,10D12.3)') (A(JJ,JJ),JJ=1,ID) CALL MATOUT (1,K,A,ID,ID) RETURN ENDIF M1=NN-KZ+1 CALL XY11 (A,I1,ID) IF(KKKK.NE.0) THEN PRINT*,'-- IN , AFTER , K,KS,I1 = ',K,KS,I1 PRINT*,'-- (INVERSE OF A, AFTER CHOLESKY) A(I,I) =' WRITE (6,'(1X,10D12.3)') (A(JJ,JJ),JJ=1,ID) CALL MATOUT (1,K,A,ID,ID) RETURN ENDIF CALL PUTAR(IADM,N,K,KS,A,ID,ID,I1,IGL) IF(M1.LT.2) GOTO 11 DO 1 IZ=2,M1 IF(IZ.GT.MB.AND.IZ.LE.M1-MR) GOTO 1 I=IZ IF(IZ.GT.M1-MR) I=MM1-1+IZ-M1 I2=IXY10(IADM,N,K,I) IF(I2.EQ.0) GOTO 1 CALL GETAR(IADM,N,K,I,B,ID,ID,IGL) CALL XY5(A,B,C,I1,I2,ID) CALL PUTAR(IADM,N,K,I,C,ID,ID,I2,IGL) 1 CONTINUE 11 I0=IXY10(IADM,N,K,MM1) IF(I0.EQ.0) GOTO 1161 CALL GETAR(IADM,N,K,MM1,B,ID,IR,IGR) CALL XY5(A,B,C,I1,I0,ID) CALL PUTAR(IADM,N,K,MM1,C,ID,IR,I0,IGR) 1161 IF(M1.LT.2) GOTO 21 M1=M1-1 IP=0 DO 22 J=1,M1 IF(J+1.GT.MB.AND.J.LE.M1-MR) GOTO 22 KJZ=KZ+J KJ=KJZ IF(KJZ.GT.NN-MR) KJ=N-NN+KJZ KS=1 IF(KJZ.GT.NN-MR) KS=MM1-1+KJZ-NN J1=J+1 IF(J1.GT.M1-MR+1) J1=MM1-1+J-M1 I3=IXY10(IADM,N,KJ,KS) I4=IXY10(IADM,N,K,J1) IF(I4.EQ.0) GOTO 22 I00=IXY10(IADM,N,KJ,MM1) IP=IP+1 CALL GETAR(IADM,N,K,J1,A,ID,ID,IGL) IF(IP.NE.1)GOTO 12 IF(I0.EQ.0) GOTO 12 IF(J .NE.1) CALL GETAR(IADM,N,K,MM1,C,ID,IR,IGR) CALL XY2(A,C,B,I1,I3,I0,ID) IF(I00.NE.0) GOTO 1162 CALL XY1 (C,I3,I0,ID,ID) GOTO 118 1162 CALL GETAR(IADM,N,KJ,MM1,C,ID,IR,IGR) 118 CALL XY6(C,B,I3,I0,ID) IF(I00.EQ.0) GOTO 1181 CALL PUTAR(IADM,N,KJ,MM1,C,ID,IR,I00,IGR) GOTO 12 1181 CALL PUTAR(IADM,N,KJ,MM1,C,ID,IR,I0,IGR) 12 N1=M1-J+1 DO 2 I=1,N1 IF(IP-1) 13,14,13 13 IIS=N1+1-I GOTO 15 14 IIS=I 15 JIS=J+IIS IF(JIS.GT.MB.AND.JIS.LE.M1+1-MR) GOTO 2 KS=1 IF(KJZ+IIS-1.GT.NN-MR) KS=MM1-2+KJZ+IIS-NN II=IIS IF(IIS.GT.NN-KJZ+1-MR) II=MM1-2+IIS-NN+KJZ JI=JIS IF(JIS.GT.M1-MR+1) JI=MM1-2+JIS-M1 KJZI=KJZ+IIS-1 IF(KJZI.GT.NN-MR) KJZI=N-NN+KJZI I6=IXY10(IADM,N,KJZI,KS) I2=IXY10(IADM,N,KJ,II) I5=IXY10(IADM,N,K,JI) IF(I5.EQ.0) GOTO 2 IF(IIS.NE.1) CALL GETAR(IADM,N,K,JI,C,ID,ID,IGL) IF(IIS-1) 18,19,18 18 CALL XY2 (A,C,B,I1,I3,I6,ID) GOTO 17 19 CALL XY8 (A,B,I1,I3,ID) 17 IF(I2.NE.0) GOTO 173 CALL XY1 (C,I3,I6,ID,ID) GOTO 174 173 CALL GETAR(IADM,N,KJ,II,C,ID,ID,IGL) 174 IF(IIS-1) 181,191,181 181 CALL XY6 (C,B,I3,I6,ID) GOTO 203 191 CALL XY7 (C,B,I3,ID) 203 CALL PUTAR(IADM,N,KJ,II,C,ID,ID,I6,IGL) 2 CONTINUE IF(IP.EQ.1)GOTO 22 IP=IP-2 IF(I0.EQ.0) GOTO 22 CALL GETAR(IADM,N,K,MM1,C,ID,IR,IGR) CALL XY2(A,C,B,I1,I3,I0,ID) IF(I00.NE.0) GOTO 204 CALL XY1 (C,I3,I0,ID,ID) GOTO 202 204 CALL GETAR(IADM,N,KJ,MM1,C,ID,IR,IGR) 202 CALL XY6(C,B,I3,I0,ID) IF(I00.EQ.0) GOTO 2021 CALL PUTAR(IADM,N,KJ,MM1,C,ID,IR,I00,IGR) GOTO 22 2021 CALL PUTAR(IADM,N,KJ,MM1,C,ID,IR,I0,IGR) 22 CONTINUE 21 CONTINUE DO 20 K=1,NN IS=NN+1-K I=IS IF(IS.GT.NN-MR) I=N-NN+IS I0=IXY10(IADM,N,I,MM1) KS=1 IF(IS.GT.NN-MR) KS=MM1-1+IS-NN I1=IXY10(IADM,N,I,KS) IF(I0.EQ.0) GOTO 251 CALL GETAR(IADM,N,I,MM1,B,ID,IR,IGR) CALL GETAR(IADM,N,I,KS,A,ID,ID,IGL) CALL XY4 (A,B,C,I1,I0,ID) 251 IF(I0.EQ.0) CALL XY1(C,I1,IR,ID,IR) CALL PUTAR(IADM,N,I,MM1,C,ID,IR,I0,IGR) M1=IS-1 IF(M1.LT.1) GOTO 20 DO 4 J=1,M1 KKZ=IS-J IF(J+1.GT.MB.AND.J.LE.NN-KKZ-MR) GOTO 4 KS=1 IF(KKZ.GT.NN-MR) KS=MM1-1+KKZ-NN KK=KKZ IF(KKZ.GT.NN-MR) KK=N-NN+KKZ J1=J+1 IF(J1.GT.NN-KKZ+1-MR) J1=MM1-1+J-NN+KKZ I2=IXY10(IADM,N,KK,KS) I3=IXY10(IADM,N,KK,J1) I03=I0*I3 I00=IXY10(IADM,N,KK,MM1) IF(I03.EQ.0) GOTO 4 CALL GETAR(IADM,N,KK,J1,A,ID,ID,IGL) CALL XY3 (A,C,B,I2,I1,I0,ID) IF(I00.NE.0) GOTO 253 CALL XY1 (A,I2,I0,ID,ID) GOTO 27 253 CALL GETAR(IADM,N,KK,MM1,A,ID,IR,IGR) 27 CALL XY6(A,B,I2,I0,ID) IF(I00.EQ.0) GOTO 271 CALL PUTAR(IADM,N,KK,MM1,A,ID,IR,I00,IGR) GOTO 4 271 CALL PUTAR(IADM,N,KK,MM1,A,ID,IR,I0,IGR) 4 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XY1(A,I1,I2,ID1,ID2) C C A=0 C DOUBLE PRECISION A,D0 DIMENSION A(ID1,ID2) DATA D0/0.0D0/ DO 1 I=1,I2 DO 1 K=1,I1 1 A(K,I)=D0 RETURN END SUBROUTINE XY2(A,B,C,L,M,N,ID) C C C=AT*B C AT=A-TRANSPOSED C DOUBLE PRECISION A,B,C DIMENSION A(1),B(1),C(1) K=1 DO 1 I=1,N CALL QTRV(A,L,M,B(K),C(K),ID) 1 K=K+ID RETURN END SUBROUTINE XY3(A,B,C,I1,I2,I3,ID) C C C=A*B C DOUBLE PRECISION A,B,C,RS,D0 DIMENSION A(ID,ID),B(ID,ID),C(ID,ID) DATA D0/0.0D0/ DO 1 I=1,I1 DO 1 K=1,I3 RS=D0 DO 2 J=1,I2 2 RS=RS+A(I,J)*B(J,K) 1 C(I,K)=RS RETURN END SUBROUTINE XY4(R,B,C,I1,I2,ID) C C C=R*B C R=UPPER TRIANGULAR MATRIX C DOUBLE PRECISION B,C,R,RS,D0 DIMENSION R(ID,ID),B(ID,ID),C(ID,ID) DATA D0/0.0D0/ DO 1 I=1,I1 DO 1 K=1,I2 RS=D0 DO 2 J=I,I1 2 RS=RS+R(I,J)*B(J,K) 1 C(I,K)=RS RETURN END SUBROUTINE XY5(R,B,C,I1,I2,ID) C C C=RT*B C RT=TRANSPOSED OF R C R=UPPER TRIANGULAR MATRIX C DOUBLE PRECISION B,C,R,RS,D0 DIMENSION R(ID,ID),B(ID,ID),C(ID,ID) DATA D0/0.0D0/ DO 1 I=1,I1 DO 1 K=1,I2 RS=D0 DO 2 J=1,I 2 RS=RS+R(J,I)*B(J,K) 1 C(I,K)=RS RETURN END SUBROUTINE XY6(A,B,I1,I2,ID) C C A=A-B C DOUBLE PRECISION A,B DIMENSION A(ID,ID),B(ID,ID) DO 1 I=1,I1 DO 1 K=1,I2 1 A(I,K)=A(I,K)-B(I,K) RETURN END SUBROUTINE XY7(A,B,I1,ID) C C A=A-B C A,B=UPPER TRIANGULAR MATRICES C DOUBLE PRECISION A,B DIMENSION A(ID,ID),B(ID,ID) DO 1 I=1,I1 DO 1 K=1,I 1 A(K,I)=A(K,I)-B(K,I) RETURN END SUBROUTINE XY8(A,R,M,N,ID) C C R=AT*A C AT=A-TRANSPOSED C ID=DIMENSIONS OF A,R C R=UPPER TRIANGULAR MATRIX C DOUBLE PRECISION A,R DIMENSION A(1),R(1) K=1 DO 1 I=1,N CALL QTRV(A,M,I,A(K),R(K),ID) 1 K=K+ID RETURN END SUBROUTINE XY9(A,N,ID) C C CHOLESKI PARTITIONING OF A C DOUBLE PRECISION A,RH,RS,D0,D1,DP40 COMMON/OUT/IOU,KKKK DIMENSION A(ID,ID) DATA D0/0.0D0/,D1/1.0D0/,DP40/1.D-40/ DO 7 K=1,N IF(K.NE.1) GOTO 2 RH=A(1,1) IF(RH.LE.DP40) CALL ERRHY IF(KKKK.NE.0) THEN PRINT*,'-- IN , ELEMENT A(1,1) = ',A(1,1) RETURN ENDIF RH=DSQRT(RH) RH=1.0D0/RH DO 1 I=1,N 1 A(1,I)=RH*A(1,I) GOTO 6 2 RH=A(K,K) IK=K-1 DO 3 J=1,IK 3 RH=-A(J,K)*A(J,K)+RH IF(RH.LE.DP40) CALL ERRHY IF(KKKK.NE.0) THEN PRINT*,'-- IN , ELEMENT ,RH,A(J,K), J,K= ',RH,A(J,K),J,K RETURN ENDIF RH=DSQRT(RH) A(K,K)=RH RH=D1/RH IK=K+1 IF(IK.GT.N) RETURN DO 5 J=IK,N RS=A(K,J) IIK=K-1 DO 4 M=1,IIK 4 RS=-A(M,J)*A(M,K)+RS 5 A(K,J)=RH*RS 6 IK=K+1 IF(IK.GT.N) RETURN DO 7 I=IK,N 7 A(I,K)=D0 RETURN END SUBROUTINE XY11(A,N,ID) C C COMPUTE INVERSE OF A C A=UPPER TRIANGULAR MATRIX C DOUBLE PRECISION A,RS,D0,D1,DP40 COMMON/OUT/IOU,KKKK DIMENSION A(ID,ID) DATA D0/0.0D0/,D1/1.0D0/,DP40/1.0D-30/ DO 1 I=1,N IF(A(I,I).LE.DP40) CALL ERRHY IF(KKKK.NE.0) THEN PRINT*,'-- IN , ELEMENT A(I,I), I= ',A(I,I),I RETURN ENDIF 1 A(I,I)=D1/A(I,I) N1=N-1 DO 2 I=1,N1 I1=I+1 DO 2 K=I1,N 2 A(K,I)=A(I,K)*A(I,I) DO 3 I=1,N1 NE=N-I DO 3 K=1,NE IK=I+K RS=D0 DO 4 L=1,I KL=K+L 4 RS=RS-A(KL,IK)*A(KL,K) 3 A(K,IK)=RS N1=N-1 DO 5 I=1,N1 I1=I+1 DO 5 K=I1,N 5 A(K,I)=D0 RETURN END SUBROUTINE QTRV(A,M,N,B,R,ID) C C R=R+B*A C DOUBLE PRECISION A,B,R,D0 DIMENSION A(ID,ID),R(ID),B(ID) DATA D0/0.0D0/ DO 1 I=1,M IF(B(I).NE.D0) GOTO 3 1 CONTINUE DO 2 I=1,N 2 R(I)=D0 RETURN 3 DO 4 K=1,N R(K)=D0 DO 4 L=I,M 4 R(K)=R(K)+B(L)*A(L,K) RETURN END SUBROUTINE ERRHY C C ERROR MESSAGES OF SUBROUTINE HYCHOL C COMMON/OUT/IOU,KKKK WRITE(IOU,1002) 1002 FORMAT(/7H >>>>> ,51HREDUCED NORMAL EQUATION MATRIX NOT POSITIVE D 1EFINIT,16X,5H<<<<<) KKKK=-1 RETURN END FUNCTION IXY10(IADM,N,I,K) C C COMPUTE RECORD NUMBER OUT OF ADDRESS IADM(I,K) C DIMENSION IADM(N,K) I2=IADM(I,K) IXY10=I2-I2/1000*1000 RETURN END SUBROUTINE GETAR(IADM,N,I,K,A,ID1,ID2,IG) C C READ RECORD WITH ADDRESS IADM(I,K) FROM FILE IG INTO A C FOR DIO-PACKAGE C DOUBLE PRECISION A DIMENSION A(ID1,ID2),IADM(N,1) I1=IADM(I,K)/1000 CALL GETR(IG,I1,A,ID1,ID2) RETURN END SUBROUTINE PUTAR(IADM,N,I,K,A,ID1,ID2,IS,IG) C C WRITE A WITH ADDRESS IADM(I,K) TO FILE IG FOR DIO-PACKAGE C DOUBLE PRECISION A COMMON/DEFIE/NOFIL,IFILE(7,2) DIMENSION A(ID1,ID2),IADM(N,1),J(1) I1=IADM(I,K)/1000 CALL SEARC(IG,IPOS) IF(IPOS.EQ.0) CALL DERRO(1,IG,J) IF(I1.NE.0) GOTO 10 I1=IFILE(7,IPOS)+1 10 IADM(I,K)=1000*I1+IS CALL PUTR(IG,I1,A,ID1,ID2) RETURN END SUBROUTINE VORHY(IADM,N,MM1,LFILE,LFNN) C C OPEN DIRECT ACCESS FILES FOR LINEAR EQUATION SYSTEM C FOR DIO-PACKAGE C CHARACTER*8 LFNN DIMENSION IADM(N,MM1),LFILE(3,2),LFNN(2) DO 1 I=1,N DO 1 K=1,MM1 1 IADM(I,K)=0 DO 2 I=1,2 2 CALL DOPEN(LFILE(1,I),LFNN(I),LFILE(2,I),LFILE(3,I),0) RETURN END SUBROUTINE EVICT(IUNIT) C C CLEAR FILE TABLE FOR FILE IUNIT FOR DIO-PACKAGE C COMMON/DEFIE/NOFIL,IFILE(7,2) DIMENSION J(1) CALL SEARC(IUNIT,IPOS) IF(IPOS.NE.0) GOTO 100 CALL DERRO(1,IUNIT,J) 100 DO 1 I=1,7 1 IFILE(I,IPOS)=0 RETURN END SUBROUTINE DERRO(NERR,IUNIT,ARRAY) C C ERROR MESSAGES OF DIO-PACKAGE C C NERR=1 UNKNOWN FILE NUMBER NERR,LFN C NERR=2 NO SPACE IN /DEFIE/ NERR,LFN C NERR=3 FILE ALREADY OPEN NERR,LFN C NERR=4 WRONG RECORD NUMBER (READ) NERR,LFN,RECNUMB C NERR=5 WRONG RECORD LENGTH (READ) NERR,LFN,RECNUMB,LRECW,LRECC C NERR=6 WRONG RECORD LENGTH (WRITE) NERR,LFN,RECNUMB,LRECW,LRECC C NERR=7 NUMBER OF RECORDS TOO SMALL NERR,LFN,RECNUMB COMMON/OUT/IOU,KKKK INTEGER ARRAY(3) IF(NERR.LT.1.OR.NERR.GT.7) GOTO 100 GOTO (1,1,1,4,5,5,7,8),NERR 1 WRITE(IOU,1000) NERR,IUNIT GOTO 200 4 WRITE(IOU,1000) NERR,IUNIT,ARRAY(1) GOTO 200 5 WRITE(IOU,1000) NERR,IUNIT,(ARRAY(I),I=1,3) GOTO 200 7 WRITE(IOU,1007) IUNIT,ARRAY(1) GOTO 200 100 WRITE(IOU,1000) NERR 200 CONTINUE STOP 8 RETURN 1000 FORMAT(/,16H >>>>> DIO-ERROR,5I10) 1007 FORMAT(/,11H >>>>> UNIT,I5,30H -- NUMBER OF RECORDS AT LEAST,I5, 123X,5H<<<<<) END SUBROUTINE SEARC(IUNIT,IPOS) C C SEARCH FILE TABLE FOR FILE IUNIT FOR DIO-PACKAGE C COMMON/DEFIE/NOFIL,IFILE(7,2) DO 1 IND=1,NOFIL IPOS=IND IF(IFILE(1,IND).EQ.IUNIT) RETURN 1 CONTINUE IPOS=0 RETURN END SUBROUTINE DINIT(MOFIL) C C INITIALIZE DIRECT ACCESS FILE TABLE FOR DIO-PACKAGE C COMMON/DEFIE/NOFIL,IFILE(7,2) DO 1 I=1,MOFIL DO 1 K=1,7 IFILE(K,I)=0 1 CONTINUE NOFIL=MOFIL RETURN C END SUBROUTINE DINIT END ------------------------------------------------MATOUT.F SUBROUTINE MATOUT (MOD,IA,A,N,DIM) C CREATION 02.02.87 C OUTPUT MATRIX. CALCULATE SPUR, MINIMAL ELEMENT.NE.0, C ALSO PRINTOUT MATRIX (1 ELEMENT = 1 CHARACTER) WHEN MOD>0 C UPDATE: LSPUR(06.03.87) IMPLICIT LOGICAL (A-Z) INTEGER MOD,IA,N,DIM REAL*8 A (DIM,DIM) INTEGER I,J,K,L,L1,L2,NN REAL*8 MIN, LSPUR, AA CHARACTER LINE*120 MIN = 9.99D66 LSPUR= 0.0D0 DO 1 I=1,N LSPUR = LSPUR + LOG ( A (I,I) ) DO 1 J=1,N AA = ABS (A(I,J)) 1 IF (AA.GT.0.0D0.AND.AA.LT.MIN) MIN = AA WRITE (6,'(1X,A,I4,2D12.3,2X,8D9.2)') + 'IA, LOG(SPUR), MIN, A(J,J) = ',IA,LSPUR,MIN, + (A(K,K),K=1,8) IF (MOD.LE.0) RETURN DO 2 NN=0, N/120 L1 = 1 + NN*120 L2 = L1 + 119 IF (L2.GT.N) L2= N PRINT*,' COLUMNS = ',L1,L2,' ', + ' 3:= >1.0D0, 2:= >1.0D-15, 1:= >1.0D-30, .:= ELSE' DO 3 I=1,DIM LINE = ' ' DO 4 L=L1,L2 J = L - NN*120 AA = ABS (A(I,L)) IF (AA.GT.1.0D0) THEN LINE (J:J) = '3' ELSEIF (AA.GT.1.0D-15) THEN LINE (J:J) = '2' ELSEIF (AA.GT.1.0D-30) THEN LINE (J:J) = '1' ELSE LINE (J:J) = '.' ENDIF 4 CONTINUE WRITE (6,'(1X,I5,2X,A)') I,LINE 3 CONTINUE 2 CONTINUE RETURN END ---------------------------------------MAIN.F PROGRAM HYC C H.B.A.P.P. LAST UPDATE 03.02.88 C C SOLUTION OF SYSTEM OF REDUCED NORMAL EQUATIONS BY SUBR. PACKAGE C H Y C H O L ((C) KLEIN). C C INITIALISE, COLLECT ALL PLATE DATA FROM -RSN- FOR INPUT IN HYCHOL C SUBROUTINES, SOLVE WHOLE REDUCED SYSTEM OF NORMAL EQUATIONS AND C OUTPUT FINAL PLATE CONSTANTS CORRECTIONS. C C DIMTI = TI AND DIMTC = TC NEEDED, EXCEPT TC = 0 THEN DIMTC = 1 C WHEN CLOSED ZONE THEN ACTUAL NUMBER OF COLUMNS IN LAST SUBMATRIX C SHOULD BE LARGER OR EQUAL TO ACTUAL NUMBER OF ALL OVERLAPPING PLATE C CONSTANTS. THIS CAN BE ACHIVED BY CHOOSING IDIM APROPRIATELY. C C +++ LOWER TRIANGEL OF REDUCED SYSTEM IS TAKEN FROM -RSN- AND C TRANSPOSED TO SUBMATRICES OF UPPER TRIANGEL FOR HYCHOL +++++ C *** NO ELEMENT BELOW DIAGONAL IS SPECIFIED FOR HYCHOL PROGRAM C PACKAGE, BECAUSE SYMMETRY IS KNOWN AND APPLIED THERE ****** C --- PGM UPGRADED FOR CLOSED ZONES (24H IN RIGHT ASCENSION) C --- IHZEI.GE.2 NEEDED C --- CHOOSE REGION OF SCRATCH FILES, FT17F001, FT18F001 BELOW C C INPUT ON 12: -RSN- REDUCED SYSTEM OF NORMAL EQ. FROM C ON 13: -B0BC- ACCUMULATED VALUES: F00, G0 C ON 5: FNRSN = FILE NAME OF -RSN- C FNB0BC = FILE NAME OF -B0BC- C FNOUTH = FILE NAME OF -OUTH- C SCALF = SCALE FACTOR TO PUSH DETERMINANT NEAR 1 C IFMOUT = .TRUE. THEN PRINTOUT OF MATRICES AB,B,C C CZ = 1 IF CLOSED ZONE (LAST IPN OVERLAP FIRST) C CMPBLA = CONNECT.MATR.PLATES, BORDER, ACTUAL LENGTH C (NUMBER OF UNKNOWNS) C TC = ACTUAL NUMBER OF COMMON PLATE CONSTANTS C IPN1,IPN2 = RANGE OF INTERNAL PLATE NUMBERS TO C PRINTOUT RESULTS (PLATE CONST. CORRECT.) C IPN1,IPN2 = ... NEW RANGE BEHIND PREVIOUS RANGE C ... (UNTIL END OF UNIT 5) C C SCRATCH ON 17: LEFT HAND SIDE OF REDUCED SYSTEM, DIRECT ACCESS C ON 18: RIGHT HAND SIDE OF REDUCED SYSTEM, DIRECT ACCESS C (17,18 ARE FILES FOR HYCHOL ONLY) C C OUTPUT ON 20: -OUTH- OUTPUT OF RESULTS FROM HYCHOL: FINAL CORRECT. C TO ALL PLATE CONSTANTS C ON 13: -B0BC- OUTPUT RESULTS FOR COMMON PLATE CONSTANTS C C PARAMETER: C DIMTI = NUMBER OF INDIVIDUAL PLATE CONSTANTS OF A PLATE C DIMTC = NUMBER OF COMMON PLATE CONSTANTS C DIMIPN = NUMBER OF PLATES IN TOTAL C MAXPC = MAXIMAL NUMBER OF PLAE CONNECTIONS C * IDIM = DIMENSION OF SUBMATIRCES USED BY HCHOL C * IHBDI = NUMBER OF SUBMATRICES IN THE DIAGONAL BAND C * IHBRA = NUMBER OF SUBMATRICES IN THE BORDER C * IHBA1 = NUMBER OF COLUMNS IN ADDRESSMATRIX = IHBDI + IHBRA + 1 C * IHZEI = NUMBER OF ROWS IN ADDRESSMATRIX C * IR = NUMBER OF RIGHT HAND SIDES (LE.IDIM REQUIRED) C * = HAVE A LOOK AT DOCUMENTS OF HYCHOL PROGRAM PACKAGE * --1.-- D E C L A R A T I O N S IMPLICIT LOGICAL (A-Z) INTEGER DIMTI,DIMTC,DIMIPN,IDIM,IHBDI,IHBRA,IHBA1,IHZEI,IR,MAXPC PARAMETER (DIMTI=8,DIMTC=1,DIMIPN=4,IDIM=16,IHBDI=2, + IHBRA=0,IHBA1=3,IHZEI=2,IR=1,MAXPC=10) * USED FOR INPUT OF -RSN-: LOGICAL LHLCZ, LHL, LHLG2, IFMOUT INTEGER I,J,K,N,II,JJ,NFL,NH,LH,IRS,IRS0,JRS,JRS0,IH0,IHL,IH,JH, + JHT,JH0,IAD,JAD,LFL,IPN,IPNC,IPNF,IDIM1,NC,CZ,TC, IFAB, + IPN0,IPN1,IPN2,INN,I0,JLAST,CMPBLA,NTOT,REST,REST0, GOMS INTEGER IPNL(MAXPC+1), IFNH(IHBDI), ILENR(IHZEI) DOUBLE PRECISION DUMMY, SCALF DOUBLE PRECISION H(IDIM,IDIM,IHBDI), FL(DIMTI,DIMTI,MAXPC+1), + F(DIMTI,DIMTI),F0(DIMTC,DIMTI),G(DIMTI),F00(DIMTC,DIMTC), + G0(DIMTC), SP(DIMTI), AB(IDIM,IDIM), AR(IDIM,IR), + A(IDIM,IDIM), B(IDIM,IDIM), C(IDIM,IDIM) CHARACTER*25 FNRSN, FNB0BC, FNOUTH * USED FOR HYCHOL SUBROUTINES: INTEGER IOU,KKKK,NOFIL,IGL,IGR,NTRSZ,LPRU, IHBDJ INTEGER IADM(IHZEI,IHBA1),LFILE(3,2) INTEGER IFILE(7,2), IXY10 * IFILE IS USED IN A COMMON BLOCK, IXY10 IS A FUNCTION CHARACTER*8 LFNAM(2) COMMON /OUT/ IOU,KKKK COMMON /IBMF/ LPRU COMMON /DEFIE/ NOFIL,IFILE * --2.-- P R E P A R A T I O N PRINT*,'- - - H. B. A. P. P. < P 8 1 H Y C > ' + ,' SOLUTION OF REDUCED NORMAL EQUATION WITH HYCHOL SUBR.' PRINT*,' ' PRINT*,'P A R A M E T E R : DIMTI, DIMTC, DIMIPN, MAXPC = ', + DIMTI, DIMTC, DIMIPN, MAXPC PRINT*,' IDIM, IHBA1, IHBDI, IHBRA, IHZEI,', + ' IR = ', IDIM, IHBA1, IHBDI, IHBRA, IHZEI, IR PRINT*,' ' READ*, FNRSN PRINT*,'INPUT FILE FNRSN = ',FNRSN READ*, FNB0BC PRINT*,'INPUT FILE FNB0BC = ',FNB0BC READ*, FNOUTH PRINT*,'OUTPUT FILE FNOUTH = ',FNOUTH OPEN (12,FILE=FNRSN ,FORM='UNFORMATTED') OPEN (13,FILE=FNB0BC,FORM='UNFORMATTED') OPEN (17,ACCESS='DIRECT',RECL=19042) OPEN (18,ACCESS='DIRECT',RECL=19042) OPEN (20,FILE=FNOUTH,FORM='UNFORMATTED') DO 21 I=1,4 21 READ (13) DUMMY READ (13) F00 READ (13) G0 PRINT*,' ' PRINT*,'PRINTOUT TO CHECK:' PRINT*,'F00 = ',F00 PRINT*,'G0 = ',G0 PRINT*,' ' READ*, SCALF PRINT*,'SCALE FACTOR SCALF = ',SCALF DO 22 I=1,DIMTC G0(I) = G0(I) * SCALF DO 22 J=1,DIMTC 22 F00(I,J) = F00(I,J) * SCALF PRINT*,' ' PRINT*,'SCALED F00 = ',F00 PRINT*,'SCALED G0 = ',G0 READ*, IFMOUT IF (IFMOUT) THEN PRINT*,'PRINTOUT OF MATRICES AB,B,C (1 ELEMENT = 1 CHARACTER)' ELSE PRINT*,'N O PRINTOUT OF MATRICES AB,B,C' ENDIF READ*, CZ IF (CZ.EQ.0) THEN PRINT*,'CZ = 0 *** NO CLOSED ZONE ***' ELSEIF (CZ.EQ.1) THEN PRINT*,'CZ = 1 *** CLOSED ZONE ***, OVERLAP FIRST/LAST PLATES' ELSE PRINT*,'CZ = ',CZ,' INVALID, USE CZ = 0 OR 1 --- S T O P ---' STOP ENDIF IF (IHZEI.LT.2) THEN PRINT*,'IHZEI IS TOO SMALL, USE IHZEI.GE.2 +++STOP+++ ' STOP ENDIF IF (IHBRA.NE.1.AND.IHBRA.NE.0) THEN PRINT*,'IHBRA = 0 OR 1 REQUIRED +++STOP+++' STOP ENDIF IHBDJ = 2 IF (IHZEI.EQ.2.AND.IHBRA.EQ.1) IHBDJ = 1 IF (IHBDI.NE.IHBDJ) THEN PRINT*,'IHBDI = ', IHBDJ,' REQUIRED +++STOP+++' STOP ENDIF PRINT*,' ' READ*, CMPBLA PRINT*,'ACTUAL LENGTH OF BORDER OF CONNECT.MAT. CMPBLA = ',CMPBLA READ*, TC PRINT*,'ACTUAL NUMBER OF COMMON PLATE CONSTANTS TC = ',TC LHL = .FALSE. LHLCZ = .FALSE. LHLG2 = .FALSE. LFL = 0 IH0 = 0 IHL = IDIM LH = 1 IDIM1 = IDIM + 1 NTOT = DIMIPN * DIMTI JLAST = NTOT - (NTOT/IDIM) * IDIM IF (JLAST.EQ.0) JLAST = IDIM PRINT*,' ' PRINT*,'ACTUAL NUMBER OF COLUMNS IN LAST SUBMATRIX (INDIVID. ', + 'PLATE CONSTANTS) JLAST = ', JLAST REST = JLAST + TC IF (CZ.NE.1) REST = TC REST0= REST - TC PRINT*,'ACTUAL NUMBER OF COLUMNS IN BORDER REST = ',REST PRINT*,'NUMBER OF UNKNOWNS IN BORDER EXEPT COMMON PLATE ', + 'CONSTANTS REST0 = ', REST0 PRINT*,' ' IF (REST.GT.IDIM) THEN PRINT*,'REST IS LARGER THAN IDIM -- INVALID -- +++STOP+++' STOP ENDIF IF (CZ.EQ.1.AND.REST0.LT.CMPBLA) THEN PRINT*,'REST IS TOO SMALL, MORE THAN ONE HYPERLINE IS NEEDED', + ' HERE TO COVER OVERLAP OF LAST PLATES WITH FIRST' PRINT*,'+++STOP+++' STOP ENDIF IFAB = 0 DO 23 NH=1,IHBDI IFNH (NH) = 0 DO 23 I=1,IDIM DO 23 J=1,IDIM H (I,J,NH) = 0.0 D0 23 CONTINUE DO 24 I = 1,IDIM DO 24 J = 1,IDIM AB (I,J) = 0.0D0 A (I,J) = 0.0D0 B (I,J) = 0.0D0 24 C (I,J) = 0.0D0 DO 25 I = 1,IDIM DO 25 J = 1,IR 25 AR (I,J) = 0.0D0 * --3.-- I N I T I A L I S E H Y C H O L * IOU= LOG.FILENUMBER FOR ERROR MESSAGES, * IGL, IGR = LOG.FILENUMBER FOR LEFT/RIGHT HAND SIDE (SCRATCH,INTERNAL) IOU = 6 IGL = 17 IGR = 18 * TRACK SIZE OF DISC IN BYTES (IBM FORMAT) NTRSZ= 19040 LPRU = NTRSZ/8 * INITIALIZATION OF FILE TABLE FOR FILES (I.E. HYCHOL ONLY 2 FILES) PRINT*,'NOW CALL DINIT' CALL DINIT(2) * INITIALIZATION FOR OPEN STATEMENTS, RIGHT AND LEFT HAND SIDE: * FILENAME, LOG.FILENUMB.,RECORD LENGTH, MAX.NUMB. OF RECORDS LFNAM (1) = 'LE.SCR' LFILE (1,1) = IGL LFILE (2,1) = IDIM*IDIM LFILE (3,1) = IHZEI*(IHBA1-1) LFNAM (2) = 'RI.SCR' LFILE (1,2) = IGR LFILE (2,2) = IDIM LFILE (3,2) = IHZEI PRINT*,'NOW CALL VORHY' CALL VORHY(IADM,IHZEI,IHBA1,LFILE,LFNAM) PRINT*,' ' * --4.-- L O O P H Y P E R L I N E S F L CALL CLOCKM (GOMS) PRINT*,'--- START LOOP INPUT CPU-TIME= ',GOMS,' MILLISEC' 4 CONTINUE DO 41 I=1,DIMTI DO 41 J=1,DIMTI DO 41 K=1,MAXPC+1 41 FL (I,J,K) = 0.0 D0 READ (12,END=900) IPN,NC,F0,G LFL = LFL + 1 DO 42 I=1,DIMTI G(I) = G(I) * SCALF DO 42 J=1,DIMTC 42 F0(J,I) = F0(J,I) * SCALF * --5.-- L O O P R E A D F MATRICES DO 5 NFL = 1,NC+1 READ (12,END=997) IPNF,IPNC,F IF (IPNF.NE.IPN) THEN PRINT*,'IPN, IPNF, IPNC, NC = ',IPN,IPNF,IPNC,NC PRINT*,'NO MATCH OF IPN, IPNF --- S T O P ---' STOP ENDIF DO 51 I=1,DIMTI DO 51 J=1,DIMTI 51 FL (I,J,NFL) = F (I,J) * SCALF IPNL (NFL) = IPNC 5 CONTINUE IRS0 = (IPN-1) * DIMTI * --6.-- L O O P L I N E S I DO 6 I = 1, DIMTI IRS = IRS0 + I IF (IRS.GT.IHL) THEN * --7.-- O U T P U T H, N E W HYPERLINE H DO 7 NH = 1, IHBDI IF (IFNH(NH).EQ.1) THEN DO 71 II=1,IDIM DO 71 JJ=1,IDIM 71 A (II,JJ) = H (JJ,II,NH) IAD = LH - NH + 1 JAD = NH IF (JAD.EQ.1) CALL MATOUT (0,IAD,A,IDIM,IDIM) CALL PUTAR (IADM,IHZEI,IAD,JAD,A,IDIM,IDIM,IDIM,IGL) ENDIF 7 CONTINUE IF (CZ.EQ.1.AND.LH.EQ.1) THEN * USE ARRAY C FOR COMMON PLATE CONSTANTS OF CLOSED ZONE * OF LINE LH = 1. WAIT FOR PLATE CONNECTIONS IN BORDER DO 72 II = 1,IDIM DO 72 JJ = REST0+1, REST0+TC 72 C(II,JJ) = AB(II,JJ) ELSEIF (CZ.EQ.1.AND.LH.EQ.IHZEI-1) THEN * USE ARRAY B FOR COMMON PLATE CONSTANTS OF CLOSED ZONE * OF LINE LH = IHZEI-1. WAIT FOR PLATE CONNECT. IN BORDER DO 73 II = 1,IDIM DO 73 JJ = REST0+1, REST0+TC 73 B(II,JJ) = AB(II,JJ) ELSEIF (REST.NE.0.AND.IFAB.EQ.1) THEN * TRANSFERE BORDER(= COMMON PLATE CONSTANTS ONLY) CALL PUTAR (IADM,IHZEI,LH ,IHBA1-1,AB,IDIM,IDIM,REST,IGL) ENDIF CALL PUTAR (IADM,IHZEI,LH ,IHBA1, AR,IDIM, IR, IR,IGR) IH0 = IHL IHL = IH0 + IDIM LH = LH + 1 IFAB = 0 IF (LH.EQ.IHZEI.AND.IHBRA.EQ.1) LHL = .TRUE. IF (LHL.AND.CZ.EQ.1) LHLCZ = .TRUE. IF (LHL.AND.IHZEI.GT.2) LHLG2 = .TRUE. DO 74 NH = 1,IHBDI IFNH (NH) = 0 DO 74 II=1,IDIM DO 74 JJ=1,IDIM 74 H (II,JJ,NH) = 0.0 D0 DO 76 II=1,IDIM DO 75 JJ=1,IDIM 75 AB (II,JJ) = 0.0 D0 DO 76 JJ=1,IR 76 AR (II,JJ) = 0.0 D0 ENDIF * - - - - - - - - (IRS.GT.IHL) -- OUTPUT HYPERLINE H - - - - - - - JH0 = IH0 + IDIM + 1 IH = IRS - IH0 * --8.-- L O O P TRANSFERE F MATRICES INTO HYPERLINE H DO 8 N = 1, NC + 1 IPNC = IPNL (N) JRS0 = (IPNC-1) * DIMTI DO 8 J=1,DIMTI JRS = JRS0 + J IF (JRS.GT.IRS) GOTO 8 JHT = JH0 - JRS NH = 1 + (JHT-1) / IDIM JH = JHT - (NH-1) * IDIM JH = IDIM1 - JH IF (LHLCZ.AND.NH.EQ.1) THEN * INTO BORDER MATRIX OF LAST HYPERLINE IFAB = 1 AB (JH,IH) = FL (I,J,N) ELSEIF (LHLG2.AND.NH.EQ.2) THEN * CONNECTIONS OF LAST PLATES WITH ADJACENT PLATES BEFORE B (JH,IH) = FL (I,J,N) ELSEIF (LHLCZ.AND.NH.EQ.IHZEI) THEN * CONNECTIONS OF LAST PLATES WITH FIRST PLATES C (JH,IH) = FL (I,J,N) ELSE IFNH (NH) = 1 H (IH,JH,NH) = FL (I,J,N) ENDIF 8 CONTINUE DO 81 J=1,TC IFAB = 1 JJ = REST0 + J 81 AB (IH,JJ) = F0 (J,I) AR (IH,1) = G (I) 6 CONTINUE * - - - - - - - - - - - - - - - - - - READ IN NEW HYPERLINE FL --- GOTO 4 * --9.-- L A S T L I N E OF HYPERMATRIX 900 CONTINUE CALL CLOCKM (GOMS) PRINT*,'--- LAST LINE OF HYPERMATRIX, CPU-TIME= ',GOMS,' MILLISEC' PRINT*,' ' PRINT*,'+ + + + LAST HYPERLINE FL READ IN + + + +' PRINT*,'LFL = ',LFL,' LH = ',LH,' IH0,IHL,IH = ',IH0,IHL,IH, + ' IPN, NC = ',IPN,NC IF (LHLCZ) THEN * LAST HYPERLINE IN CASE OF CLOSED ZONE * INPUT F00, G0 INTO BORDER AB AND RIGHT HAND SIDE AR: DO 91 I=1,TC II = REST0 + I AR (II,1) = G0 (I) DO 91 J=1,TC JJ = REST0 + J 91 AB (II,JJ) = F00 (J,I) IF (IFMOUT) THEN PRINT*,'MATOUT: ARRAY AB: AT POSITION = ',IHZEI,IHBA1-1 CALL MATOUT (1,IHZEI,AB,JLAST,IDIM) PRINT*,'MATOUT: ARRAY C : AT POSITION = 1 ',IHBA1-1 CALL MATOUT (1, 1,C,JLAST,IDIM) PRINT*,'MATOUT: ARRAY B : AT POSITION = ',IHZEI-1,IHBA1-1 CALL MATOUT (1,IHZEI-1,B,JLAST,IDIM) ENDIF CALL PUTAR (IADM,IHZEI,IHZEI ,IHBA1 ,AR,IDIM, IR, IR,IGR) CALL PUTAR (IADM,IHZEI,IHZEI ,IHBA1-1,AB,IDIM,IDIM,REST,IGL) CALL PUTAR (IADM,IHZEI, 1,IHBA1-1,C ,IDIM,IDIM,REST,IGL) IF (IHZEI.GT.2) + CALL PUTAR (IADM,IHZEI,IHZEI-1,IHBA1-1,B ,IDIM,IDIM,REST,IGL) ELSE * LAST HYPERLINE IN CASE OF NO CLOSED ZONE IF (IH.GT.0) THEN DO 92 NH = 1, IHBDI IF (IFNH(NH).EQ.1) THEN * SUBMATRIX A IS NOT EMTY, TRANSFERE TO HYCHOL DO 93 II=1,IDIM DO 93 JJ=1,IDIM 93 A (II,JJ) = H (JJ,II,NH) IAD = LH - NH + 1 JAD = NH CALL PUTAR (IADM,IHZEI,IAD,JAD,A,IDIM,IDIM,JLAST,IGL) ENDIF 92 CONTINUE ENDIF IF (REST.NE.0) THEN DO 94 I=1,TC II = REST0 + I AR (II,1) = G0 (I) DO 94 J=1,TC JJ = REST0 + J 94 AB (II,JJ) = F00 (J,I) CALL PUTAR (IADM,IHZEI,IHZEI,IHBA1-1,AB,IDIM,IDIM,REST,IGL) CALL PUTAR (IADM,IHZEI,IHZEI,IHBA1, AR,IDIM, IR, IR,IGR) ENDIF ENDIF * --10.-- S T A R T S U B R O U T I N E H C H O L CALL CLOCKM (GOMS) PRINT*,'--- START H C H O L CPU-TIME= ',GOMS,' MILLISEC' * FORMATION OF LINEAR EQUATION SYSTEM FAILED IF... IF(KKKK.NE.0) GOTO 1000 PRINT*,' ' PRINT*,'*******************************************************' PRINT*,'INPUT FROM - R S N - FINISHED START H C H O L' PRINT*,'*******************************************************' PRINT*,' ' CALL HCHOL (IADM,IHZEI,IHBA1,IHZEI,IHBDI,IHBRA,A,B,C, + IDIM,IR,IGL,IGR) CALL CLOCKM (GOMS) PRINT*,'--- AFTER H C H O L CPU-TIME= ',GOMS,' MILLISEC' * SOLUTION OF LINEAR EQUATION SYSTEM FAILED IF... IF(KKKK.NE.0) GOTO 1000 PRINT*,'SOLUTION OF LINEAR EQUTION SYSTEM FINISHED + + + + + + +' * --11.-- G E T S O L U T I O N IPN = 0 I = 0 INN = IHZEI - IHBRA I0 = IHZEI - IHBA1 + 1 DO 111 N = 1,DIMTI 111 SP(N) = 0.0 D0 DO 11 N = 1, IHZEI II = 0 CALL GETAR (IADM,IHZEI,N,IHBA1,A,IDIM,IR,IGR) * GET ACTUAL LENGTH OF SOLUTION VECTOR OF ROW N K = 1 IF (N.GT.INN) K = N - I0 ILENR (N) = IXY10 (IADM,IHZEI,N,K) 110 CONTINUE II = II + 1 IF (II.GT.ILENR(N)) GOTO 11 I = I + 1 IF (I.GT.DIMTI) THEN * ALL PLATE CONSTANTS OF A PLATE ARE COLLECTED, OUTPUT TO -OUTH- IPN = IPN + 1 WRITE (20) IPN, SP DO 112 JJ=1,DIMTI 112 SP (JJ) = 0.0 D0 I = 1 ENDIF SP (I) = A (II,1) GOTO 110 11 CONTINUE IF (TC.EQ.0) THEN PRINT*,'NO COMMON PLATE CONSTATNTS, LAST PLATE CONSTANTS = ' IPN = IPN + 1 WRITE (6,222) IPN,SP ELSE IPN = 0 PRINT*,'LAST PLATE CONSTANTS = COMMON PLATE CONSTANTS :' WRITE (6,222) IPN,SP 222 FORMAT (1X,I5,9F14.10) WRITE (13) IPN,SP ENDIF WRITE (20) IPN, SP CLOSE (20) * --12.-- T E S T P R I N T O U T CALL CLOCKM (GOMS) PRINT*,'--- AFTER SOLUTION CPU-TIME= ',GOMS,' MILLISEC' PRINT*,' ' PRINT*,'********** T E S T P R I N T O U T ***********' OPEN (20,FILE=FNOUTH,FORM='UNFORMATTED') IPN0 = 0 12 CONTINUE READ (5,*,END=999) IPN1,IPN2 PRINT*,' ' PRINT*,'N E W R A N G E : IPN1, IPN2 = ',IPN1,IPN2 PRINT*,'============================================' DO 121 I = IPN0+1, IPN1-1 121 READ (20,END=998) J DO 122 I = IPN1, IPN2 READ (20,END=998) IPN, SP WRITE (6,222) IPN, SP 122 CONTINUE IPN0 = IPN2 GOTO 12 997 PRINT*,' ' PRINT*,'END OF FILE ON 12 DURING READ OF NC+1 F MATRICES' PRINT*,'LFL, NFL, NC = ',LFL,NFL,NC GOTO 999 998 PRINT*,' ' PRINT*,'END OF FILE ON 20 DURING READ OF -OUTH- , IPN = ',I GOTO 999 1000 PRINT*,' ' PRINT*,' E R R O R * * * GOTO 1000 S T O P, KKKK = ',KKKK 999 STOP END SUBROUTINE MATOUT (MOD,IA,A,N,DIM) C CREATION 02.02.87 C OUTPUT MATRIX. CALCULATE SPUR, MINIMAL ELEMENT.NE.0, C ALSO PRINTOUT MATRIX (1 ELEMENT = 1 CHARACTER) WHEN MOD>0 C UPDATE: LSPUR(06.03.87) IMPLICIT LOGICAL (A-Z) INTEGER MOD,IA,N,DIM REAL*8 A (DIM,DIM) INTEGER I,J,K,L,L1,L2,NN REAL*8 MIN, LSPUR, AA CHARACTER LINE*120 MIN = 9.99D66 LSPUR= 0.0D0 DO 1 I=1,N LSPUR = LSPUR + LOG ( A (I,I) ) DO 1 J=1,N AA = ABS (A(I,J)) 1 IF (AA.GT.0.0D0.AND.AA.LT.MIN) MIN = AA WRITE (6,'(1X,A,I4,2D12.3,2X,8D9.2)') + 'IA, LOG(SPUR), MIN, A(J,J) = ',IA,LSPUR,MIN, + (A(K,K),K=1,8) IF (MOD.LE.0) RETURN DO 2 NN=0, N/120 L1 = 1 + NN*120 L2 = L1 + 119 IF (L2.GT.N) L2= N PRINT*,' COLUMNS = ',L1,L2,' ', + ' 3:= >1.0D0, 2:= >1.0D-15, 1:= >1.0D-30, .:= ELSE' DO 3 I=1,DIM LINE = ' ' DO 4 L=L1,L2 J = L - NN*120 AA = ABS (A(I,L)) IF (AA.GT.1.0D0) THEN LINE (J:J) = '3' ELSEIF (AA.GT.1.0D-15) THEN LINE (J:J) = '2' ELSEIF (AA.GT.1.0D-30) THEN LINE (J:J) = '1' ELSE LINE (J:J) = '.' ENDIF 4 CONTINUE WRITE (6,'(1X,I5,2X,A)') I,LINE 3 CONTINUE 2 CONTINUE RETURN END SUBROUTINE CLOCKM(GOMS) C SHOULD RETURN TIME IN SECONDS, BUT WE DON'T KNOW HOW TO DO IT (CAH) INTEGER GOMS GOMS=0 RETURN END -----------------------------------------------end message