diff --git a/docker/Dockerfile b/docker/Dockerfile
new file mode 100644
index 0000000000000000000000000000000000000000..d3daa1c5fd94e6453951c7bdd71c2c6acb14e180
--- /dev/null
+++ b/docker/Dockerfile
@@ -0,0 +1,110 @@
+FROM continuumio/anaconda:latest
+MAINTAINER Pierre-Elouan Rethore <pe@retho.re>
+
+RUN apt-get update \
+ && apt-get install -y \
+ build-essential \
+ gcc \
+ gfortran \
+ wget \
+ binutils \
+ gcc \
+ libgtk2.0-0 \
+ libgtk2.0-dev \
+ psmisc
+
+
+
+RUN mkdir /install
+WORKDIR /install
+
+#ENV REQUIREMENTS requirements1.txt
+#ARG REQUIREMENTS
+
+# Install python library requirements to run the notebooks
+COPY requirements2.txt /install
+RUN pip install --upgrade pip \
+ && pip install -r /install/requirements2.txt
+
+## Install PyOptSparse & IpOpt
+RUN apt-get update \
+ && apt-get install -y \
+ mercurial \
+ meld
+RUN apt-get update \
+ && apt-get install -y \
+ swig
+
+ENV POSDIR /install/pyoptsparse
+ENV IPV 3.11.7
+ENV IPOPT_DIR $POSDIR/pyoptsparse/pyIPOPT/Ipopt
+
+
+RUN hg clone https://bitbucket.org/mdolab/pyoptsparse $POSDIR
+
+# Install Ipopt
+COPY install_ipopt2.sh /install
+COPY ma27ad.f /install
+RUN /install/install_ipopt2.sh
+ENV LD_LIBRARY_PATH $LD_LIBRARY_PATH:$IPOPT_DIR/lib
+
+# Install SNOPT
+COPY snopt/* $POSDIR/pyoptsparse/pySNOPT/source/
+
+## Install PyOptSparse
+#COPY install_pyoptsparse.sh /install
+WORKDIR $POSDIR
+RUN python setup.py install
+
+#RUN mkdir /notebooks
+#WORKDIR /notebooks
+
+#COPY mycert.pem /install/
+#COPY mykey.key /install/
+#COPY jupyter_notebook_config.py /install/
+
+# Add Tini. Tini operates as a process subreaper for jupyter. This prevents
+# kernel crashes.
+#ENV TINI_VERSION v0.6.0
+#ADD https://github.com/krallin/tini/releases/download/${TINI_VERSION}/tini /usr/bin/tini
+#RUN chmod +x /usr/bin/tini
+
+
+# Install the Colonel
+
+RUN mkdir /deb
+WORKDIR /deb
+COPY *.deb /deb/
+RUN dpkg -i *.deb
+
+RUN apt-get clean \
+ && apt-get autoremove -y
+
+RUN mkdir /install
+WORKDIR /install
+COPY fuga/*.pas /install/
+COPY fuga/*.lpr /install/
+COPY fuga/*.lpi /install/
+
+
+## Build
+#RUN lazbuild ColonelLazarus.lpr
+
+#RUN curl -sL https://deb.nodesource.com/setup_8.x | bash -
+#RUN apt-get update -y && apt-get install -y nodejs
+
+
+#RUN jupyter labextension install @jupyter-widgets/jupyterlab-manager \
+# && jupyter labextension install jupyterlab_bokeh
+
+#ENTRYPOINT ["/usr/bin/tini", "--"]
+
+
+#EXPOSE 8898
+
+CMD bash
+
+ # CMD jupyter notebook \
+ # --notebook-dir=/notebooks \
+ # --config=/install/jupyter_notebook_config.py \
+ # --allow-root
diff --git a/docker/fpc-src_3.0.4-2_amd64.deb b/docker/fpc-src_3.0.4-2_amd64.deb
new file mode 100644
index 0000000000000000000000000000000000000000..0c121d0fa64ace642134c8559cce2c341eb9e6a8
Binary files /dev/null and b/docker/fpc-src_3.0.4-2_amd64.deb differ
diff --git a/docker/fpc_3.0.4-2_amd64.deb b/docker/fpc_3.0.4-2_amd64.deb
new file mode 100644
index 0000000000000000000000000000000000000000..986469225241e34ee0228b4d60c5d6eef6cbe4c1
Binary files /dev/null and b/docker/fpc_3.0.4-2_amd64.deb differ
diff --git a/docker/install_ipopt2.sh b/docker/install_ipopt2.sh
new file mode 100755
index 0000000000000000000000000000000000000000..6ae360eb659aa76acde920f229863cc4ef6eb027
--- /dev/null
+++ b/docker/install_ipopt2.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+curl https://www.coin-or.org/download/source/Ipopt/Ipopt-$IPV.tgz > /install/Ipopt-$IPV.tgz
+cd /install
+tar xvf Ipopt-$IPV.tgz
+mv Ipopt-$IPV $IPOPT_DIR
+cp /install/ma27ad.f $IPOPT_DIR/ThirdParty/HSLold/
+cd $IPOPT_DIR/ThirdParty/Blas
+sh ./get.Blas
+cd $IPOPT_DIR/ThirdParty/Lapack
+sh ./get.Lapack
+cd $IPOPT_DIR
+./configure --disable-linear-solver-loader
+make install
diff --git a/docker/ma27ad.f b/docker/ma27ad.f
new file mode 100644
index 0000000000000000000000000000000000000000..d5ea17a9f3757eea3052dc806463b884ffc230aa
--- /dev/null
+++ b/docker/ma27ad.f
@@ -0,0 +1,3496 @@
+* COPYRIGHT (c) 1982 AEA Technology
+*######DATE 20 September 2001
+C September 2001: threadsafe version of MA27
+C 19/3/03. Array ICNTL in MA27GD made assumed size.
+C 17/9/09. Bug corrected in MA27HD. Also some tidying of code in
+C MA27HD and change of most comments to lower case.
+C 16/6/10. Statement function in MA27OD replaced by in-line code.
+
+ SUBROUTINE MA27ID(ICNTL,CNTL)
+ INTEGER ICNTL(30)
+ DOUBLE PRECISION CNTL(5)
+
+ INTEGER IFRLVL
+ PARAMETER ( IFRLVL=5 )
+
+C Stream number for error messages
+ ICNTL(1) = 6
+C Stream number for diagnostic messages
+ ICNTL(2) = 6
+C Control the level of diagnostic printing.
+C 0 no printing
+C 1 printing of scalar parameters and first parts of arrays.
+C 2 printing of scalar parameters and whole of arrays.
+ ICNTL(3) = 0
+C The largest integer such that all integers I in the range
+C -ICNTL(4).LE.I.LE.ICNTL(4) can be handled by the shortest integer
+C type in use.
+ ICNTL(4) = 2139062143
+C Minimum number of eliminations in a step that is automatically
+C accepted. if two adjacent steps can be combined and each has less
+C eliminations then they are combined.
+ ICNTL(5) = 1
+C Control whether direct or indirect access is used by MA27C/CD.
+C Indirect access is employed in forward and back substitution
+C respectively if the size of a block is less than
+C ICNTL(IFRLVL+MIN(10,NPIV)) and ICNTL(IFRLVL+10+MIN(10,NPIV))
+C respectively, where NPIV is the number of pivots in the block.
+ ICNTL(IFRLVL+1) = 32639
+ ICNTL(IFRLVL+2) = 32639
+ ICNTL(IFRLVL+3) = 32639
+ ICNTL(IFRLVL+4) = 32639
+ ICNTL(IFRLVL+5) = 14
+ ICNTL(IFRLVL+6) = 9
+ ICNTL(IFRLVL+7) = 8
+ ICNTL(IFRLVL+8) = 8
+ ICNTL(IFRLVL+9) = 9
+ ICNTL(IFRLVL+10) = 10
+ ICNTL(IFRLVL+11) = 32639
+ ICNTL(IFRLVL+12) = 32639
+ ICNTL(IFRLVL+13) = 32639
+ ICNTL(IFRLVL+14) = 32689
+ ICNTL(IFRLVL+15) = 24
+ ICNTL(IFRLVL+16) = 11
+ ICNTL(IFRLVL+17) = 9
+ ICNTL(IFRLVL+18) = 8
+ ICNTL(IFRLVL+19) = 9
+ ICNTL(IFRLVL+20) = 10
+C Not used
+ ICNTL(26) = 0
+ ICNTL(27) = 0
+ ICNTL(28) = 0
+ ICNTL(29) = 0
+ ICNTL(30) = 0
+
+C Control threshold pivoting.
+ CNTL(1) = 0.1D0
+C If a column of the reduced matrix has relative density greater than
+C CNTL(2), it is forced into the root. All such columns are taken to
+C have sparsity pattern equal to their merged patterns, so the fill
+C and operation counts may be overestimated.
+ CNTL(2) = 1.0D0
+C An entry with absolute value less than CNTL(3) is never accepted as
+C a 1x1 pivot or as the off-diagonal of a 2x2 pivot.
+ CNTL(3) = 0.0D0
+C Not used
+ CNTL(4) = 0.0
+ CNTL(5) = 0.0
+
+ RETURN
+ END
+
+ SUBROUTINE MA27AD(N,NZ,IRN,ICN,IW,LIW,IKEEP,IW1,NSTEPS,IFLAG,
+ + ICNTL,CNTL,INFO,OPS)
+C THIS SUBROUTINE COMPUTES A MINIMUM DEGREE ORDERING OR ACCEPTS A GIVEN
+C ORDERING. IT COMPUTES ASSOCIATED ASSEMBLY AND ELIMINATION
+C INFORMATION FOR MA27B/BD.
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT
+C ALTERED.
+C IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE
+C NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
+C TO IW (SEE DESCRIPTION OF IW).
+C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE
+C NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
+C TO IW (SEE DESCRIPTION OF IW).
+C IW NEED NOT BE SET ON INPUT. IT IS USED FOR WORKSPACE.
+C IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE
+C EQUIVALENCED TO IW(K), WHERE K.GT.NZ.
+C LIW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+3*N
+C FOR THE IFLAG=0 ENTRY AND AT LEAST NZ+3*N FOR THE IFLAG=1
+C ENTRY. IT IS NOT ALTERED.
+C IKEEP NEED NOT BE SET UNLESS AN ORDERING IS GIVEN, IN WHICH CASE
+C IKEEP(I,1) MUST BE SET TO THE POSITION OF VARIABLE I IN THE
+C ORDER. ON OUTPUT IKEEP CONTAINS INFORMATION NEEDED BY MA27B/BD.
+C IKEEP(I,1) HOLDS THE POSITION OF VARIABLE I IN THE PIVOT ORDER.
+C IKEEP(I,2), IKEEP(I,3) HOLD THE NUMBER OF ELIMINATIONS, ASSEMBLIES
+C AT MAJOR STEP I, I=1,2,...,NSTEPS. NOTE THAT WHEN AN ORDER IS
+C GIVEN IT MAY BE REPLACED BY ANOTHER ORDER THAT GIVES IDENTICAL
+C NUMERICAL RESULTS.
+C IW1 IS USED FOR WORKSPACE.
+C NSTEPS NEED NOT BE SET. ON OUTPUT IT CONTAINS THE NUMBER OF MAJOR
+C STEPS NEEDED FOR A DEFINITE MATRIX AND MUST BE PASSED UNCHANGED
+C TO MA27B/BD.
+C IFLAG MUST SET TO ZERO IF THE USER WANTS THE PIVOT ORDER CHOSEN
+C AUTOMATICALLY AND TO ONE IF HE WANTS TO SPECIFY IT IN IKEEP.
+C ICNTL is an INTEGER array of length 30 containing user options
+C with integer values (defaults set in MA27I/ID)
+C ICNTL(1) (LP) MUST BE SET TO THE STREAM NUMBER FOR ERROR MESSAGES.
+C ERROR MESSAGES ARE SUPPRESSED IF ICNTL(1) IS NOT POSITIVE.
+C IT IS NOT ALTERED.
+C ICNTL(2) (MP) MUST BE SET TO THE STREAM NUMBER FOR DIAGNOSTIC
+C MESSAGES. DIAGNOSTIC MESSAGES ARE SUPPRESSED IF ICNTL(2) IS NOT
+C POSITIVE. IT IS NOT ALTERED.
+C ICNTL(3) (LDIAG) CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING.
+C 0 NO PRINTING
+C 1 PRINTING OF SCALAR PARAMETERS AND FIRST PARTS OF ARRAYS.
+C 2 PRINTING OF SCALAR PARAMETERS AND WHOLE OF ARRAYS.
+C ICNTL(4) (IOVFLO) IS THE LARGEST INTEGER SUCH THAT ALL INTEGERS
+C I IN THE RANGE -IOVFLO.LE.I.LE.IOVFLO CAN BE HANDLED BY THE
+C SHORTEST INTEGER TYPE IN USE.
+C ICNT(5) (NEMIN) MUST BE SET TO THE MINIMUM NUMBER OF ELIMINATIONS
+C IN A STEP THAT IS AUTOMATICALLY ACCEPTED. IF TWO ADJACENT STEPS
+C CAN BE COMBINED AND EACH HAS LESS ELIMINATIONS THEN THEY ARE
+C COMBINED.
+C ICNTL(IFRLVL+I) I=1,20, (IFRLVL) MUST BE SET TO CONTROL WHETHER
+C DIRECT OR INDIRECT ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS
+C IS EMPLOYED IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE
+C SIZE OF A BLOCK IS LESS THAN ICNTL(IFRLVL+(MIN(10,NPIV)) AND
+C ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
+C NUMBER OF PIVOTS IN THE BLOCK.
+C ICNTL(I) I=26,30 are not used.
+C CNTL is an DOUBLE PRECISION array of length 5 containing user options
+C with real values (defaults set in MA27I/ID)
+C CNTL(1) (U) IS USED TO CONTROL THRESHOLD PIVOTING. IT IS NOT
+C ALTERED.
+C CNTL(2) (FRATIO) has default value 1.0. If a column of the
+C reduced matrix has relative density greater than CNTL(2), it
+C is forced into the root. All such columns are taken to have
+C sparsity pattern equal to their merged patterns, so the fill
+C and operation counts may be overestimated.
+C CNTL(3) (PIVTOL) has default value 0.0. An entry with absolute
+C value less than CNTL(3) is never accepted as a 1x1 pivot or
+C as the off-diagonal of a 2x2 pivot.
+C CNTL(I) I=4,5 are not used.
+C INFO is an INTEGER array of length 20 which is used to return
+C information to the user.
+C INFO(1) (IFLAG) is an error return code, zero for success, greater
+C than zero for a warning and less than zero for errors, see
+C INFO(2).
+C INFO(2) (IERROR) HOLDS ADDITIONAL INFORMATION IN THE EVENT OF ERRORS.
+C IF INFO(1)=-3 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR IW.
+C IF INFO(1)=-4 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR A.
+C IF INFO(1)=-5 INFO(2) IS SET TO THE PIVOT STEP AT WHICH SINGULARITY
+C WAS DETECTED.
+C IF INFO(1)=-6 INFO(2) IS SET TO THE PIVOT STEP AT WHICH A CHANGE OF
+C PIVOT SIGN WAS FOUND.
+C IF INFO(1)= 1 INFO(2) HOLDS THE NUMBER OF FAULTY ENTRIES.
+C IF INFO(1)= 2 INFO(2) IS SET TO THE NUMBER OF SIGNS CHANGES IN
+C THE PIVOTS.
+C IF INFO(1)=3 INFO(2) IS SET TO THE RANK OF THE MATRIX.
+C INFO(3) and INFO(4) (NRLTOT and NIRTOT) REAL AND INTEGER STRORAGE
+C RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF NO COMPRESSES ARE
+C ALLOWED.
+C INFO(5) and INFO(6) (NRLNEC and NIRNEC) REAL AND INTEGER STORAGE
+C RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF COMPRESSES ARE
+C ALLOWED AND THE MATRIX IS DEFINITE.
+C INFO(7) and INFO(8) (NRLADU and NIRADU) REAL AND INTEGER STORAGE
+C RESPECTIVELY FOR THE MATRIX FACTORS AS CALCULATED BY MA27A/AD
+C FOR THE DEFINITE CASE.
+C INFO(9) and INFO(10) (NRLBDU and NIRBDU) REAL AND INTEGER STORAGE
+C RESPECTIVELY FOR THE MATRIX FACTORS AS FOUND BY MA27B/BD.
+C INFO(11) (NCMPA) ACCUMULATES THE NUMBER OF TIMES THE ARRAY IW IS
+C COMPRESSED BY MA27A/AD.
+C INFO(12) and INFO(13) (NCMPBR and NCMPBI) ACCUMULATE THE NUMBER
+C OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY MA27B/BD.
+C INFO(14) (NTWO) IS USED BY MA27B/BD TO RECORD THE NUMBER OF 2*2
+C PIVOTS USED.
+C INFO(15) (NEIG) IS USED BY ME27B/BD TO RECORD THE NUMBER OF
+C NEGATIVE EIGENVALUES OF A.
+C INFO(16) to INFO(20) are not used.
+C OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE
+C TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE.
+C
+C .. Scalar Arguments ..
+ INTEGER IFLAG,LIW,N,NSTEPS,NZ
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION CNTL(5),OPS
+ INTEGER ICNTL(30),INFO(20)
+ INTEGER ICN(*),IKEEP(N,3),IRN(*),IW(LIW),IW1(N,2)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IWFR,K,L1,L2,LLIW
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27GD,MA27HD,MA27JD,MA27KD,MA27LD,MA27MD,MA27UD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MIN
+C ..
+C .. Executable Statements ..
+ DO 5 I = 1,15
+ INFO(I) = 0
+ 5 CONTINUE
+
+ IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 40
+C PRINT INPUT VARIABLES.
+ WRITE (ICNTL(2),FMT=10) N,NZ,LIW,IFLAG
+
+ 10 FORMAT(/,/,' ENTERING MA27AD WITH N NZ LIW IFLAG',
+ + /,21X,I7,I7,I9,I7)
+
+ NSTEPS = 0
+ K = MIN(8,NZ)
+ IF (ICNTL(3).GT.1) K = NZ
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=20) (IRN(I),ICN(I),I=1,K)
+
+ 20 FORMAT (' MATRIX NON-ZEROS',/,4 (I9,I6),/,
+ + (I9,I6,I9,I6,I9,I6,I9,I6))
+
+ K = MIN(10,N)
+ IF (ICNTL(3).GT.1) K = N
+ IF (IFLAG.EQ.1 .AND. K.GT.0) THEN
+ WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)
+ END IF
+
+ 30 FORMAT (' IKEEP(.,1)=',10I6,/, (12X,10I6))
+
+ 40 IF (N.LT.1 .OR. N.GT.ICNTL(4)) GO TO 70
+ IF (NZ.LT.0) GO TO 100
+ LLIW = LIW - 2*N
+ L1 = LLIW + 1
+ L2 = L1 + N
+ IF (IFLAG.EQ.1) GO TO 50
+ IF (LIW.LT.2*NZ+3*N+1) GO TO 130
+C SORT
+ CALL MA27GD(N,NZ,IRN,ICN,IW,LLIW,IW1,IW1(1,2),IW(L1),IWFR,
+ + ICNTL,INFO)
+C ANALYZE USING MINIMUM DEGREE ORDERING
+ CALL MA27HD(N,IW1,IW,LLIW,IWFR,IW(L1),IW(L2),IKEEP(1,2),
+ + IKEEP(1,3),IKEEP,ICNTL(4),INFO(11),CNTL(2))
+ GO TO 60
+C SORT USING GIVEN ORDER
+ 50 IF (LIW.LT.NZ+3*N+1) GO TO 120
+ CALL MA27JD(N,NZ,IRN,ICN,IKEEP,IW,LLIW,IW1,IW1(1,2),IW(L1),IWFR,
+ + ICNTL,INFO)
+C ANALYZE USING GIVEN ORDER
+ CALL MA27KD(N,IW1,IW,LLIW,IWFR,IKEEP,IKEEP(1,2),IW(L1),IW(L2),
+ + INFO(11))
+C PERFORM DEPTH-FIRST SEARCH OF ASSEMBLY TREE
+ 60 CALL MA27LD(N,IW1,IW(L1),IKEEP,IKEEP(1,2),IKEEP(1,3),IW(L2),
+ + NSTEPS,ICNTL(5))
+C EVALUATE STORAGE AND OPERATION COUNT REQUIRED BY MA27B/BD IN THE
+C DEFINITE CASE.
+C SET IW(1) SO THAT ARRAYS IW AND IRN CAN BE TESTED FOR EQUIVALENCE
+C IN MA27M/MD.
+ IF(NZ.GE.1) IW(1) = IRN(1) + 1
+ CALL MA27MD(N,NZ,IRN,ICN,IKEEP,IKEEP(1,3),IKEEP(1,2),IW(L2),
+ + NSTEPS,IW1,IW1(1,2),IW,INFO,OPS)
+ GO TO 160
+
+ 70 INFO(1) = -1
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+
+ 80 FORMAT (' **** ERROR RETURN FROM MA27AD **** INFO(1)=',I3)
+
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=90) N
+
+ 90 FORMAT (' VALUE OF N OUT OF RANGE ... =',I10)
+
+ GO TO 160
+
+ 100 INFO(1) = -2
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=110) NZ
+
+ 110 FORMAT (' VALUE OF NZ OUT OF RANGE .. =',I10)
+
+ GO TO 160
+
+ 120 INFO(2) = NZ + 3*N + 1
+ GO TO 140
+
+ 130 INFO(2) = 2*NZ + 3*N + 1
+ 140 INFO(1) = -3
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=150) LIW,INFO(2)
+
+ 150 FORMAT (' LIW TOO SMALL, MUST BE INCREASED FROM',I10,
+ + ' TO AT LEAST',I10)
+
+ 160 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0 .OR. INFO(1).LT.0) GO TO 200
+C PRINT PARAMETER VALUES ON EXIT.
+ WRITE (ICNTL(2),FMT=170) NSTEPS,INFO(1),OPS,INFO(2),INFO(3),
+ + INFO(4),INFO(5),INFO(6),INFO(7),INFO(8),INFO(11)
+
+ 170 FORMAT (/,' LEAVING MA27AD WITH NSTEPS INFO(1) OPS IERROR',
+ + ' NRLTOT NIRTOT',
+ + /,20X,2I7,F7.0,3I7,
+ + /,20X,' NRLNEC NIRNEC NRLADU NIRADU NCMPA',
+ + /,20X,6I7)
+
+ K = MIN(9,N)
+ IF (ICNTL(3).GT.1) K = N
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)
+ K = MIN(K,NSTEPS)
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=180) (IKEEP(I,2),I=1,K)
+
+ 180 FORMAT (' IKEEP(.,2)=',10I6,/, (12X,10I6))
+
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=190) (IKEEP(I,3),I=1,K)
+
+ 190 FORMAT (' IKEEP(.,3)=',10I6,/, (12X,10I6))
+
+ 200 CONTINUE
+
+ RETURN
+ END
+
+ SUBROUTINE MA27BD(N,NZ,IRN,ICN,A,LA,IW,LIW,IKEEP,NSTEPS,MAXFRT,
+ + IW1,ICNTL,CNTL,INFO)
+C THIS SUBROUTINE COMPUTES THE FACTORISATION OF THE MATRIX INPUT IN
+C A,IRN,ICN USING INFORMATION (IN IKEEP) FROM MA27A/AD.
+C N MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED.
+C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT
+C ALTERED.
+C IRN,ICN,A. ENTRY K (K=1,NZ) OF IRN,ICN MUST BE SET TO THE ROW
+C AND COLUMN INDEX RESPECTIVELY OF THE NON-ZERO IN A.
+C IRN AND ICN ARE UNALTERED BY MA27B/BD.
+C ON EXIT, ENTRIES 1 TO NRLBDU OF A HOLD REAL INFORMATION
+C ON THE FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD.
+C LA LENGTH OF ARRAY A. AN INDICATION OF A SUITABLE VALUE,
+C SUFFICIENT FOR FACTORIZATION OF A DEFINITE MATRIX, WILL
+C HAVE BEEN PROVIDED IN NRLNEC AND NRLTOT BY MA27A/AD.
+C IT IS NOT ALTERED BY MA27B/BD.
+C IW NEED NOT BE SET ON ENTRY. USED AS A WORK ARRAY BY MA27B/BD.
+C ON EXIT, ENTRIES 1 TO NIRBDU HOLD INTEGER INFORMATION ON THE
+C FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD.
+C LIW LENGTH OF ARRAY IW. AN INDICATION OF A SUITABLE VALUE WILL
+C HAVE BEEN PROVIDED IN NIRNEC AND NIRTOT BY MA27A/AD.
+C IT IS NOT ALTERED BY MA27B/BD.
+C IKEEP MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
+C IT IS NOT ALTERED BY MA27B/BD.
+C NSTEPS MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
+C IT IS NOT ALTERED BY MA27B/BD.
+C MAXFRT NEED NOT BE SET AND MUST BE PASSED UNCHANGED TO MA27C/CD.
+C IT IS THE MAXIMUM SIZE OF THE FRONTAL MATRIX GENERATED DURING
+C THE DECOMPOSITION.
+C IW1 USED AS WORKSPACE BY MA27B/BD.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C CNTL is a DOUBLE PRECISION array of length 5, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER LA,LIW,MAXFRT,N,NSTEPS,NZ
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA),CNTL(5)
+ INTEGER ICN(*),IKEEP(N,3),IRN(*),IW(LIW),IW1(N)
+ INTEGER ICNTL(30),INFO(20)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IAPOS,IBLK,IPOS,IROWS,J1,J2,JJ,K,KBLK,KZ,LEN,NCOLS,
+ + NROWS,NZ1
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27ND,MA27OD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MIN
+C ..
+C .. Executable Statements ..
+ INFO(1) = 0
+
+ IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 60
+C PRINT INPUT PARAMETERS.
+ WRITE (ICNTL(2),FMT=10) N,NZ,LA,LIW,NSTEPS,CNTL(1)
+
+ 10 FORMAT (/,/,
+ + ' ENTERING MA27BD WITH N NZ LA LIW',
+ + ' NSTEPS U',/,21X,I7,I7,I9,I9,I7,1PD10.2)
+
+ KZ = MIN(6,NZ)
+ IF (ICNTL(3).GT.1) KZ = NZ
+ IF (NZ.GT.0) WRITE (ICNTL(2),FMT=20) (A(K),IRN(K),ICN(K),K=1,KZ)
+
+ 20 FORMAT (' MATRIX NON-ZEROS',/,1X,2 (1P,D16.3,2I6),/,
+ + (1X,1P,D16.3,2I6,1P,D16.3,2I6))
+
+ K = MIN(9,N)
+ IF (ICNTL(3).GT.1) K = N
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)
+
+ 30 FORMAT (' IKEEP(.,1)=',10I6,/, (12X,10I6))
+
+ K = MIN(K,NSTEPS)
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=40) (IKEEP(I,2),I=1,K)
+
+ 40 FORMAT (' IKEEP(.,2)=',10I6,/, (12X,10I6))
+
+ IF (K.GT.0) WRITE (ICNTL(2),FMT=50) (IKEEP(I,3),I=1,K)
+
+ 50 FORMAT (' IKEEP(.,3)=',10I6,/, (12X,10I6))
+
+ 60 IF (N.LT.1 .OR. N.GT.ICNTL(4)) GO TO 70
+ IF (NZ.LT.0) GO TO 100
+ IF (LIW.LT.NZ) GO TO 120
+ IF (LA.LT.NZ+N) GO TO 150
+ IF (NSTEPS.LT.1 .OR. NSTEPS.GT.N) GO TO 175
+C SORT
+ CALL MA27ND(N,NZ,NZ1,A,LA,IRN,ICN,IW,LIW,IKEEP,IW1,ICNTL,INFO)
+ IF (INFO(1).EQ.-3) GO TO 130
+ IF (INFO(1).EQ.-4) GO TO 160
+C FACTORIZE
+ CALL MA27OD(N,NZ1,A,LA,IW,LIW,IKEEP,IKEEP(1,3),NSTEPS,MAXFRT,
+ + IKEEP(1,2),IW1,ICNTL,CNTL,INFO)
+ IF (INFO(1).EQ.-3) GO TO 130
+ IF (INFO(1).EQ.-4) GO TO 160
+ IF (INFO(1).EQ.-5) GO TO 180
+ IF (INFO(1).EQ.-6) GO TO 200
+C **** WARNING MESSAGE ****
+ IF (INFO(1).EQ.3 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=65) INFO(1),INFO(2)
+ END IF
+
+ 65 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
+ + ' *** INFO(1) =',I2,
+ + /,5X,'MATRIX IS SINGULAR. RANK=',I5)
+
+ GO TO 220
+C **** ERROR RETURNS ****
+ 70 INFO(1) = -1
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+
+ 80 FORMAT (' **** ERROR RETURN FROM MA27BD **** INFO(1)=',I3)
+
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=90) N
+
+ 90 FORMAT (' VALUE OF N OUT OF RANGE ... =',I10)
+
+ GO TO 220
+
+ 100 INFO(1) = -2
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=110) NZ
+
+ 110 FORMAT (' VALUE OF NZ OUT OF RANGE .. =',I10)
+
+ GO TO 220
+
+ 120 INFO(1) = -3
+ INFO(2) = NZ
+ 130 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=140) LIW,INFO(2)
+
+ 140 FORMAT (' LIW TOO SMALL, MUST BE INCREASED FROM',I10,' TO',
+ + ' AT LEAST',I10)
+
+ GO TO 220
+
+ 150 INFO(1) = -4
+ INFO(2) = NZ + N
+ 160 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=170) LA,INFO(2)
+
+ 170 FORMAT (' LA TOO SMALL, MUST BE INCREASED FROM ',I10,' TO',
+ + ' AT LEAST',I10)
+
+ GO TO 220
+
+ 175 INFO(1) = -7
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) THEN
+ WRITE (ICNTL(1),FMT='(A)') ' NSTEPS is out of range'
+ END IF
+ GO TO 220
+
+ 180 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=190) INFO(2)
+
+ 190 FORMAT (' ZERO PIVOT AT STAGE',I10,
+ + ' WHEN INPUT MATRIX DECLARED DEFINITE')
+
+ GO TO 220
+
+ 200 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
+ IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=210)
+
+ 210 FORMAT (' CHANGE IN SIGN OF PIVOT ENCOUNTERED',
+ + ' WHEN FACTORING ALLEGEDLY DEFINITE MATRIX')
+
+ 220 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0 .OR. INFO(1).LT.0) GO TO 310
+C PRINT OUTPUT PARAMETERS.
+ WRITE (ICNTL(2),FMT=230) MAXFRT,INFO(1),INFO(9),INFO(10),INFO(12),
+ + INFO(13),INFO(14),INFO(2)
+
+ 230 FORMAT (/,' LEAVING MA27BD WITH',
+ + /,10X,' MAXFRT INFO(1) NRLBDU NIRBDU NCMPBR',
+ + ' NCMPBI NTWO IERROR',
+ + /,11X,8I7)
+
+ IF (INFO(1).LT.0) GO TO 310
+C PRINT OUT MATRIX FACTORS FROM MA27B/BD.
+ KBLK = ABS(IW(1)+0)
+ IF (KBLK.EQ.0) GO TO 310
+ IF (ICNTL(3).EQ.1) KBLK = 1
+ IPOS = 2
+ IAPOS = 1
+ DO 300 IBLK = 1,KBLK
+ NCOLS = IW(IPOS)
+ NROWS = IW(IPOS+1)
+ J1 = IPOS + 2
+ IF (NCOLS.GT.0) GO TO 240
+ NCOLS = -NCOLS
+ NROWS = 1
+ J1 = J1 - 1
+ 240 WRITE (ICNTL(2),FMT=250) IBLK,NROWS,NCOLS
+
+ 250 FORMAT (' BLOCK PIVOT =',I8,' NROWS =',I8,' NCOLS =',I8)
+
+ J2 = J1 + NCOLS - 1
+ IPOS = J2 + 1
+ WRITE (ICNTL(2),FMT=260) (IW(JJ),JJ=J1,J2)
+
+ 260 FORMAT (' COLUMN INDICES =',10I6,/, (17X,10I6))
+
+ WRITE (ICNTL(2),FMT=270)
+
+ 270 FORMAT (' REAL ENTRIES .. EACH ROW STARTS ON A NEW LINE')
+
+ LEN = NCOLS
+ DO 290 IROWS = 1,NROWS
+ J1 = IAPOS
+ J2 = IAPOS + LEN - 1
+ WRITE (ICNTL(2),FMT=280) (A(JJ),JJ=J1,J2)
+
+ 280 FORMAT (1P,5D13.3)
+
+ LEN = LEN - 1
+ IAPOS = J2 + 1
+ 290 CONTINUE
+ 300 CONTINUE
+ 310 RETURN
+ END
+
+ SUBROUTINE MA27CD(N,A,LA,IW,LIW,W,MAXFRT,RHS,IW1,NSTEPS,
+ + ICNTL,INFO)
+C THIS SUBROUTINE USES THE FACTORISATION OF THE MATRIX IN A,IW TO
+C SOLVE A SYSTEM OF EQUATIONS.
+C N MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED.
+C A,IW HOLD INFORMATION ON THE FACTORS AND MUST BE UNCHANGED SINCE
+C THE CALL TO MA27B/BD. THEY ARE NOT ALTERED BY MA27C/CDD.
+C LA,LIW MUST BE SET TO THE LENGTHS OF A,IW RESPECTIVELY. THEY
+C ARE NOT ALTERED.
+C W USED AS A WORK ARRAY.
+C MAXFRT IS THE LENGTH OF W AND MUST BE PASSED UNCHANGED FROM THE
+C CALL TO MA27B/BD. IT IS NOT ALTERED.
+C RHS MUST BE SET TO THE RIGHT HAND SIDE FOR THE EQUATIONS BEING
+C SOLVED. ON EXIT, THIS ARRAY WILL HOLD THE SOLUTION.
+C IW1 USED AS A WORK ARRAY.
+C NSTEPS IS THE LENGTH OF IW1 WHICH MUST BE AT LEAST THE ABSOLUTE
+C VALUE OF IW(1). IT IS NOT ALTERED.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER LA,LIW,MAXFRT,N,NSTEPS
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA),RHS(N),W(MAXFRT)
+ INTEGER IW(LIW),IW1(NSTEPS),ICNTL(30),INFO(20)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IAPOS,IBLK,IPOS,IROWS,J1,J2,JJ,K,KBLK,LATOP,LEN,NBLK,
+ + NCOLS,NROWS
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27QD,MA27RD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MIN
+C ..
+C .. Executable Statements ..
+ INFO(1) = 0
+
+ IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 110
+C PRINT INPUT PARAMETERS.
+ WRITE (ICNTL(2),FMT=10) N,LA,LIW,MAXFRT,NSTEPS
+
+ 10 FORMAT (/,/,' ENTERING MA27CD WITH N LA LIW MAXFRT',
+ + ' NSTEPS',/,21X,5I7)
+C PRINT OUT MATRIX FACTORS FROM MA27B/BD.
+ KBLK = ABS(IW(1)+0)
+ IF (KBLK.EQ.0) GO TO 90
+ IF (ICNTL(3).EQ.1) KBLK = 1
+ IPOS = 2
+ IAPOS = 1
+ DO 80 IBLK = 1,KBLK
+ NCOLS = IW(IPOS)
+ NROWS = IW(IPOS+1)
+ J1 = IPOS + 2
+ IF (NCOLS.GT.0) GO TO 20
+ NCOLS = -NCOLS
+ NROWS = 1
+ J1 = J1 - 1
+ 20 WRITE (ICNTL(2),FMT=30) IBLK,NROWS,NCOLS
+
+ 30 FORMAT (' BLOCK PIVOT =',I8,' NROWS =',I8,' NCOLS =',I8)
+
+ J2 = J1 + NCOLS - 1
+ IPOS = J2 + 1
+ WRITE (ICNTL(2),FMT=40) (IW(JJ),JJ=J1,J2)
+
+ 40 FORMAT (' COLUMN INDICES =',10I6,/, (17X,10I6))
+
+ WRITE (ICNTL(2),FMT=50)
+
+ 50 FORMAT (' REAL ENTRIES .. EACH ROW STARTS ON A NEW LINE')
+
+ LEN = NCOLS
+ DO 70 IROWS = 1,NROWS
+ J1 = IAPOS
+ J2 = IAPOS + LEN - 1
+ WRITE (ICNTL(2),FMT=60) (A(JJ),JJ=J1,J2)
+
+ 60 FORMAT (1P,5D13.3)
+
+ LEN = LEN - 1
+ IAPOS = J2 + 1
+ 70 CONTINUE
+ 80 CONTINUE
+ 90 K = MIN(10,N)
+ IF (ICNTL(3).GT.1) K = N
+ IF (N.GT.0) WRITE (ICNTL(2),FMT=100) (RHS(I),I=1,K)
+
+ 100 FORMAT (' RHS',1P,5D13.3,/, (4X,1P,5D13.3))
+
+ 110 IF (IW(1).LT.0) GO TO 130
+ NBLK = IW(1)
+ IF (NBLK.GT.0) GO TO 140
+C WE HAVE A ZERO MATRIX
+ DO 120 I = 1,N
+ RHS(I) = 0.0D0
+ 120 CONTINUE
+ GO TO 150
+
+ 130 NBLK = -IW(1)
+C FORWARD SUBSTITUTION
+ 140 CALL MA27QD(N,A,LA,IW(2),LIW-1,W,MAXFRT,RHS,IW1,NBLK,LATOP,ICNTL)
+C BACK SUBSTITUTION.
+ CALL MA27RD(N,A,LA,IW(2),LIW-1,W,MAXFRT,RHS,IW1,NBLK,LATOP,ICNTL)
+ 150 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 170
+C PRINT OUTPUT PARAMETERS.
+ WRITE (ICNTL(2),FMT=160)
+
+ 160 FORMAT (/,/,' LEAVING MA27CD WITH')
+
+ IF (N.GT.0) WRITE (ICNTL(2),FMT=100) (RHS(I),I=1,K)
+ 170 CONTINUE
+
+ RETURN
+ END
+ SUBROUTINE MA27GD(N,NZ,IRN,ICN,IW,LW,IPE,IQ,FLAG,IWFR,
+ + ICNTL,INFO)
+C
+C SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27H/HD.
+C
+C GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC
+C MATRIX, CONSTRUCT THE SPARSITY PATTERN OF THE OFF-DIAGONAL
+C PART OF THE WHOLE MATRIX (UPPER AND LOWER TRIANGULAR PARTS).
+C EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT
+C THE PAIR. DIAGONAL ELEMENTS AND DUPLICATES ARE IGNORED.
+C
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT
+C ALTERED.
+C IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW NUMBERS OF THE
+C NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
+C TO IW (SEE DESCRIPTION OF IW).
+C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE
+C NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
+C TO IW (SEE DESCRIPTION OF IW).
+C IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF
+C COLUMN INDICES, EACH LIST BEING HEADED BY ITS LENGTH.
+C IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE
+C EQUIVALENCED TO IW(K), WHERE K.GT.NZ.
+C LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+N.
+C IT IS NOT ALTERED.
+C IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF
+C THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY.
+C IQ NEED NOT BE SET. ON OUTPUT IQ(I),I=1,N CONTAINS THE NUMBER OF
+C OFF-DIAGONAL N0N-ZEROS IN ROW I INCLUDING DUPLICATES.
+C FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE ENTRIES
+C TO BE IDENTIFIED QUICKLY.
+C IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST
+C UNUSED LOCATION IN IW.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER IWFR,LW,N,NZ
+C ..
+C .. Array Arguments ..
+ INTEGER FLAG(N),ICN(*),IPE(N),IQ(N),IRN(*),IW(LW)
+ INTEGER ICNTL(*),INFO(20)
+C ..
+C .. Local Scalars ..
+ INTEGER I,ID,J,JN,K,K1,K2,L,LAST,LR,N1,NDUP
+C ..
+C .. Executable Statements ..
+C
+C INITIALIZE INFO(2) AND COUNT IN IPE THE
+C NUMBERS OF NON-ZEROS IN THE ROWS AND MOVE ROW AND COLUMN
+C NUMBERS INTO IW.
+ INFO(2) = 0
+ DO 10 I = 1,N
+ IPE(I) = 0
+ 10 CONTINUE
+ LR = NZ
+ IF (NZ.EQ.0) GO TO 120
+ DO 110 K = 1,NZ
+ I = IRN(K)
+ J = ICN(K)
+ IF (I.LT.J) THEN
+ IF (I.GE.1 .AND. J.LE.N) GO TO 90
+ ELSE IF (I.GT.J) THEN
+ IF (J.GE.1 .AND. I.LE.N) GO TO 90
+ ELSE
+ IF (I.GE.1 .AND. I.LE.N) GO TO 80
+ END IF
+ INFO(2) = INFO(2) + 1
+ INFO(1) = 1
+ IF (INFO(2).LE.1 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=60) INFO(1)
+ END IF
+
+ 60 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27AD',
+ + ' *** INFO(1) =',I2)
+
+ IF (INFO(2).LE.10 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=70) K,I,J
+ END IF
+
+ 70 FORMAT (I6,'TH NON-ZERO (IN ROW',I6,' AND COLUMN',I6,
+ + ') IGNORED')
+
+ 80 I = 0
+ J = 0
+ GO TO 100
+
+ 90 IPE(I) = IPE(I) + 1
+ IPE(J) = IPE(J) + 1
+ 100 IW(K) = J
+ LR = LR + 1
+ IW(LR) = I
+ 110 CONTINUE
+C
+C ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW STARTS IN BOTH IPE AND IQ
+C AND INITIALIZE FLAG
+ 120 IQ(1) = 1
+ N1 = N - 1
+ IF (N1.LE.0) GO TO 140
+ DO 130 I = 1,N1
+ FLAG(I) = 0
+ IF (IPE(I).EQ.0) IPE(I) = -1
+ IQ(I+1) = IPE(I) + IQ(I) + 1
+ IPE(I) = IQ(I)
+ 130 CONTINUE
+ 140 LAST = IPE(N) + IQ(N)
+ FLAG(N) = 0
+ IF (LR.GE.LAST) GO TO 160
+ K1 = LR + 1
+ DO 150 K = K1,LAST
+ IW(K) = 0
+ 150 CONTINUE
+ 160 IPE(N) = IQ(N)
+ IWFR = LAST + 1
+ IF (NZ.EQ.0) GO TO 230
+C
+C RUN THROUGH PUTTING THE MATRIX ELEMENTS IN THE RIGHT PLACE
+C BUT WITH SIGNS INVERTED. IQ IS USED FOR HOLDING RUNNING POINTERS
+C AND IS LEFT HOLDING POINTERS TO ROW ENDS.
+ DO 220 K = 1,NZ
+ J = IW(K)
+ IF (J.LE.0) GO TO 220
+ L = K
+ IW(K) = 0
+ DO 210 ID = 1,NZ
+ IF (L.GT.NZ) GO TO 170
+ L = L + NZ
+ GO TO 180
+
+ 170 L = L - NZ
+ 180 I = IW(L)
+ IW(L) = 0
+ IF (I.LT.J) GO TO 190
+ L = IQ(J) + 1
+ IQ(J) = L
+ JN = IW(L)
+ IW(L) = -I
+ GO TO 200
+
+ 190 L = IQ(I) + 1
+ IQ(I) = L
+ JN = IW(L)
+ IW(L) = -J
+ 200 J = JN
+ IF (J.LE.0) GO TO 220
+ 210 CONTINUE
+ 220 CONTINUE
+C
+C RUN THROUGH RESTORING SIGNS, REMOVING DUPLICATES AND SETTING THE
+C MATE OF EACH NON-ZERO.
+C NDUP COUNTS THE NUMBER OF DUPLICATE ELEMENTS.
+ 230 NDUP = 0
+ DO 280 I = 1,N
+ K1 = IPE(I) + 1
+ K2 = IQ(I)
+ IF (K1.LE.K2) GO TO 240
+C ROW IS EMPTY. SET POINTER TO ZERO.
+ IPE(I) = 0
+ IQ(I) = 0
+ GO TO 280
+C ON ENTRY TO THIS LOOP FLAG(J).LT.I FOR J=1,2,...,N. DURING THE LOOP
+C FLAG(J) IS SET TO I IF A NON-ZERO IN COLUMN J IS FOUND. THIS
+C PERMITS DUPLICATES TO BE RECOGNIZED QUICKLY.
+ 240 DO 260 K = K1,K2
+ J = -IW(K)
+ IF (J.LE.0) GO TO 270
+ L = IQ(J) + 1
+ IQ(J) = L
+ IW(L) = I
+ IW(K) = J
+ IF (FLAG(J).NE.I) GO TO 250
+ NDUP = NDUP + 1
+ IW(L) = 0
+ IW(K) = 0
+ 250 FLAG(J) = I
+ 260 CONTINUE
+ 270 IQ(I) = IQ(I) - IPE(I)
+ IF (NDUP.EQ.0) IW(K1-1) = IQ(I)
+ 280 CONTINUE
+ IF (NDUP.EQ.0) GO TO 310
+C
+C COMPRESS IW TO REMOVE DUMMY ENTRIES CAUSED BY DUPLICATES.
+ IWFR = 1
+ DO 300 I = 1,N
+ K1 = IPE(I) + 1
+ IF (K1.EQ.1) GO TO 300
+ K2 = IQ(I) + IPE(I)
+ L = IWFR
+ IPE(I) = IWFR
+ IWFR = IWFR + 1
+ DO 290 K = K1,K2
+ IF (IW(K).EQ.0) GO TO 290
+ IW(IWFR) = IW(K)
+ IWFR = IWFR + 1
+ 290 CONTINUE
+ IW(L) = IWFR - L - 1
+ 300 CONTINUE
+ 310 RETURN
+
+ END
+
+ SUBROUTINE MA27HD(N,IPE,IW,LW,IWFR,NV,NXT,LST,IPD,FLAG,IOVFLO,
+ + NCMPA,FRATIO)
+C Bug was found in the then identical subroutine MA57HD which was then
+C corrected.
+C Changes made in September 2009 because of bug in compress control
+C found by Nick.
+C
+C ANALYSIS SUBROUTINE
+C
+C GIVEN REPRESENTATION OF THE WHOLE MATRIX (EXCLUDING DIAGONAL)
+C PERFORM MINIMUM DEGREE ORDERING, CONSTRUCTING TREE POINTERS.
+C IT WORKS WITH SUPERVARIABLES WHICH ARE COLLECTIONS OF ONE OR MORE
+C VARIABLES, STARTING WITH SUPERVARIABLE I CONTAINING VARIABLE I FOR
+C I = 1,2,...,N. ALL VARIABLES IN A SUPERVARIABLE ARE ELIMINATED
+C TOGETHER. EACH SUPERVARIABLE HAS AS NUMERICAL NAME THAT OF ONE
+C OF ITS VARIABLES (ITS PRINCIPAL VARIABLE).
+C
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE
+C START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF-
+C DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS. IF
+C SUPERVARIABLE I IS ABSORBED INTO SUPERVARIABLE J THEN IPE(I)=-J.
+C IF SUPERVARIABLE I IS ELIMINATED THEN IPE(I) EITHER POINTS TO THE
+C LIST OF SUPERVARIABLES FOR CREATED ELEMENT I OR IS ZERO IF
+C THE CREATED ELEMENT IS NULL. IF ELEMENT I
+C IS ABSORBED INTO ELEMENT J THEN IPE(I)=-J.
+C IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY
+C ROWS, EACH LIST BEING HEADED BY ITS LENGTH.
+C DURING EXECUTION THESE LISTS ARE REVISED AND HOLD
+C LISTS OF ELEMENTS AND SUPERVARIABLES. THE ELEMENTS
+C ALWAYS HEAD THE LISTS. WHEN A SUPERVARIABLE
+C IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF SUPERVARIABLES
+C IN THE NEW ELEMENT.
+C LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED.
+C IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE.
+C IT IS REVISED DURING EXECUTION AND CONTINUES TO HAVE THIS MEANING.
+C NV(JS) NEED NOT BE SET. DURING EXECUTION IT IS ZERO IF
+C JS IS NOT A PRINCIPAL VARIABLE AND IF IT IS IT HOLDS
+C THE NUMBER OF VARIABLES IN SUPERVARIABLE JS. FOR ELIMINATED
+C VARIABLES IT IS SET TO THE DEGREE AT THE TIME OF ELIMINATION.
+C NXT(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE NEXT
+C SUPERVARIABLE HAVING THE SAME DEGREE AS JS, OR ZERO
+C IF IT IS LAST IN ITS LIST.
+C LST(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE
+C LAST SUPERVARIABLE HAVING THE SAME DEGREE AS JS OR
+C -(ITS DEGREE) IF IT IS FIRST IN ITS LIST.
+C IPD(ID) NEED NOT BE SET. DURING EXECUTION IT
+C IS THE FIRST SUPERVARIABLE WITH DEGREE ID OR ZERO
+C IF THERE ARE NONE.
+C FLAG IS USED AS WORKSPACE FOR ELEMENT AND SUPERVARIABLE FLAGS.
+C WHILE THE CODE IS FINDING THE DEGREE OF SUPERVARIABLE IS
+C FLAG HAS THE FOLLOWING VALUES.
+C A) FOR THE CURRENT PIVOT/NEW ELEMENT ME
+C FLAG(ME)=-1
+C B) FOR VARIABLES JS
+C FLAG(JS)=-1 IF JS IS NOT A PRINCIPAL VARIABLE
+C FLAG(JS)=0 IF JS IS A SUPERVARIABLE IN THE NEW ELEMENT
+C FLAG(JS)=NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW
+C ELEMENT THAT HAS BEEN COUNTED IN THE DEGREE
+C CALCULATION
+C FLAG(JS).GT.NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW
+C ELEMENT THAT HAS NOT BEEN COUNTED IN THE DEGREE
+C CALCULATION
+C C) FOR ELEMENTS IE
+C FLAG(IE)=-1 IF ELEMENT IE HAS BEEN MERGED INTO ANOTHER
+C FLAG(IE)=-NFLG IF ELEMENT IE HAS BEEN USED IN THE DEGREE
+C CALCULATION FOR IS.
+C FLAG(IE).LT.-NFLG IF ELEMENT IE HAS NOT BEEN USED IN THE
+C DEGREE CALCULATION FOR IS
+C IOVFLO see ICNTL(4) in MA27A/AD.
+C NCMPA see INFO(11) in MA27A/AD.
+C FRATIO see CNTL(2) in MA27A/AD.
+C
+C .. Scalar Arguments ..
+ DOUBLE PRECISION FRATIO
+ INTEGER IWFR,LW,N,IOVFLO,NCMPA
+C ..
+C .. Array Arguments ..
+ INTEGER FLAG(N),IPD(N),IPE(N),IW(LW),LST(N),NV(N),NXT(N)
+C ..
+C .. Local Scalars ..
+ INTEGER I,ID,IDL,IDN,IE,IP,IS,JP,JP1,JP2,JS,K,K1,K2,KE,KP,KP0,KP1,
+ + KP2,KS,L,LEN,LIMIT,LN,LS,LWFR,MD,ME,ML,MS,NEL,NFLG,NP,
+ + NP0,NS,NVPIV,NVROOT,ROOT
+C LIMIT Limit on number of variables for putting node in root.
+C NVROOT Number of variables in the root node
+C ROOT Index of the root node (N+1 if none chosen yet).
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27UD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MIN
+C ..
+C If a column of the reduced matrix has relative density greater than
+C CNTL(2), it is forced into the root. All such columns are taken to
+C have sparsity pattern equal to their merged patterns, so the fill
+C and operation counts may be overestimated.
+C
+C IS,JS,KS,LS,MS,NS are used to refer to supervariables.
+C IE,JE,KE are used to refer to elements.
+C IP,JP,KP,K,NP are used to point to lists of elements
+C or supervariables.
+C ID is used for the degree of a supervariable.
+C MD is used for the current minimum degree.
+C IDN is used for the no. of variables in a newly created element
+C NEL is used to hold the no. of variables that have been
+C eliminated.
+C ME=MS is the name of the supervariable eliminated and
+C of the element created in the main loop.
+C NFLG is used for the current flag value in array FLAG. It starts
+C with the value IOVFLO and is reduced by 1 each time it is used
+C until it has the value 2 when it is reset to the value IOVFLO.
+C
+C .. Executable Statements ..
+C Initializations
+ DO 10 I = 1,N
+ IPD(I) = 0
+ NV(I) = 1
+ FLAG(I) = IOVFLO
+ 10 CONTINUE
+ MD = 1
+ NCMPA = 0
+ NFLG = IOVFLO
+ NEL = 0
+ ROOT = N+1
+ NVROOT = 0
+C
+C Link together variables having same degree
+ DO 30 IS = 1,N
+ K = IPE(IS)
+ IF (K.GT.0) THEN
+ ID = IW(K) + 1
+ NS = IPD(ID)
+ IF (NS.GT.0) LST(NS) = IS
+ NXT(IS) = NS
+ IPD(ID) = IS
+ LST(IS) = -ID
+ ELSE
+C We have a variable that can be eliminated at once because there is
+C no off-diagonal nonzero in its row.
+ NEL = NEL + 1
+ FLAG(IS) = -1
+ NXT(IS) = 0
+ LST(IS) = 0
+ ENDIF
+ 30 CONTINUE
+
+C
+C Start of main loop
+C
+ DO 340 ML = 1,N
+C Leave loop if all variables have been eliminated.
+ IF (NEL+NVROOT+1.GE.N) GO TO 350
+C
+C Find next supervariable for elimination.
+ DO 40 ID = MD,N
+ MS = IPD(ID)
+ IF (MS.GT.0) GO TO 50
+ 40 CONTINUE
+ 50 MD = ID
+C Nvpiv holds the number of variables in the pivot.
+ NVPIV = NV(MS)
+C
+C Remove chosen variable from linked list
+ NS = NXT(MS)
+ NXT(MS) = 0
+ LST(MS) = 0
+ IF (NS.GT.0) LST(NS) = -ID
+ IPD(ID) = NS
+ ME = MS
+ NEL = NEL + NVPIV
+C IDN holds the degree of the new element.
+ IDN = 0
+C
+C Run through the list of the pivotal supervariable, setting tree
+C pointers and constructing new list of supervariables.
+C KP is a pointer to the current position in the old list.
+ KP = IPE(ME)
+ FLAG(MS) = -1
+C IP points to the start of the new list.
+ IP = IWFR
+C LEN holds the length of the list associated with the pivot.
+ LEN = IW(KP)
+ DO 140 KP1 = 1,LEN
+ KP = KP + 1
+ KE = IW(KP)
+C Jump if KE is an element that has not been merged into another.
+ IF (FLAG(KE).LE.-2) GO TO 60
+C Jump if KE is an element that has been merged into another or is
+C a supervariable that has been eliminated.
+ IF (FLAG(KE).LE.0) THEN
+ IF (IPE(KE).NE.-ROOT) GO TO 140
+C KE has been merged into the root
+ KE = ROOT
+ IF (FLAG(KE).LE.0) GO TO 140
+ END IF
+C We have a supervariable. Prepare to search rest of list.
+ JP = KP - 1
+ LN = LEN - KP1 + 1
+ IE = MS
+ GO TO 70
+C Search variable list of element KE, using JP as a pointer to it.
+ 60 IE = KE
+ JP = IPE(IE)
+ LN = IW(JP)
+C
+C Search for different supervariables and add them to the new list,
+C compressing when necessary. This loop is executed once for
+C each element in the list and once for all the supervariables
+C in the list.
+ 70 DO 130 JP1 = 1,LN
+ JP = JP + 1
+ IS = IW(JP)
+C Jump if IS is not a principal variable or has already been counted.
+ IF (FLAG(IS).LE.0) THEN
+ IF (IPE(IS).EQ.-ROOT) THEN
+C IS has been merged into the root
+ IS = ROOT
+ IW(JP) = ROOT
+ IF (FLAG(IS).LE.0) GO TO 130
+ ELSE
+ GO TO 130
+ END IF
+ END IF
+ FLAG(IS) = 0
+
+C To fix Nick bug need to add one here to store (eventually) length
+C of new row
+ IF (IWFR .GE. LW-1) THEN
+C Logic was previously as below
+CCC IF (IWFR.LT.LW) GO TO 100
+C Prepare for compressing IW by adjusting pointers and
+C lengths so that the lists being searched in the inner and outer
+C loops contain only the remaining entries.
+ IPE(MS) = KP
+ IW(KP) = LEN - KP1
+ IPE(IE) = JP
+ IW(JP) = LN - JP1
+C Compress IW
+ CALL MA27UD(N,IPE,IW,IP-1,LWFR,NCMPA)
+C Copy new list forward
+ JP2 = IWFR - 1
+ IWFR = LWFR
+ IF (IP.GT.JP2) GO TO 90
+ DO 80 JP = IP,JP2
+ IW(IWFR) = IW(JP)
+ IWFR = IWFR + 1
+ 80 CONTINUE
+C Adjust pointers for the new list and the lists being searched.
+ 90 IP = LWFR
+ JP = IPE(IE)
+ KP = IPE(ME)
+ ENDIF
+
+C Store IS in new list.
+ IW(IWFR) = IS
+ IDN = IDN + NV(IS)
+ IWFR = IWFR + 1
+C Remove IS from degree linked list
+ LS = LST(IS)
+ LST(IS) = 0
+ NS = NXT(IS)
+ NXT(IS) = 0
+ IF (NS.GT.0) LST(NS) = LS
+ IF (LS.LT.0) THEN
+ LS = -LS
+ IPD(LS) = NS
+ ELSE IF (LS.GT.0) THEN
+ NXT(LS) = NS
+ END IF
+ 130 CONTINUE
+C Jump if we have just been searching the variables at the end of
+C the list of the pivot.
+ IF (IE.EQ.MS) GO TO 150
+C Set tree pointer and flag to indicate element IE is absorbed into
+C new element ME.
+ IPE(IE) = -ME
+ FLAG(IE) = -1
+ 140 CONTINUE
+
+C Store the degree of the pivot.
+ 150 NV(MS) = IDN + NVPIV
+
+C Jump if new element is null.
+ IF (IWFR.EQ.IP) THEN
+ IPE(ME) = 0
+ GO TO 340
+ ENDIF
+
+ K1 = IP
+ K2 = IWFR - 1
+C
+C Run through new list of supervariables revising each associated list,
+C recalculating degrees and removing duplicates.
+ LIMIT = NINT(FRATIO*(N-NEL))
+ DO 310 K = K1,K2
+ IS = IW(K)
+ IF (IS.EQ.ROOT) GO TO 310
+ IF (NFLG.GT.2) GO TO 170
+C Reset FLAG values to +/-IOVFLO.
+ DO 160 I = 1,N
+ IF (FLAG(I).GT.0) FLAG(I) = IOVFLO
+ IF (FLAG(I).LE.-2) FLAG(I) = -IOVFLO
+ 160 CONTINUE
+ NFLG = IOVFLO
+C Reduce NFLG by one to cater for this supervariable.
+ 170 NFLG = NFLG - 1
+C Begin with the degree of the new element. Its variables must always
+C be counted during the degree calculation and they are already
+C flagged with the value 0.
+ ID = IDN
+C Run through the list associated with supervariable IS
+ KP1 = IPE(IS) + 1
+C NP points to the next entry in the revised list.
+ NP = KP1
+ KP2 = IW(KP1-1) + KP1 - 1
+ DO 220 KP = KP1,KP2
+ KE = IW(KP)
+C Test whether KE is an element, a redundant entry or a supervariable.
+ IF (FLAG(KE).EQ.-1) THEN
+ IF (IPE(KE).NE.-ROOT) GO TO 220
+C KE has been merged into the root
+ KE = ROOT
+ IW(KP) = ROOT
+ IF (FLAG(KE).EQ.-1) GO TO 220
+ END IF
+ IF (FLAG(KE).GE.0) GO TO 230
+C Search list of element KE, revising the degree when new variables
+C found.
+ JP1 = IPE(KE) + 1
+ JP2 = IW(JP1-1) + JP1 - 1
+ IDL = ID
+ DO 190 JP = JP1,JP2
+ JS = IW(JP)
+C Jump if JS has already been counted.
+ IF (FLAG(JS).LE.NFLG) GO TO 190
+ ID = ID + NV(JS)
+ FLAG(JS) = NFLG
+ 190 CONTINUE
+C Jump if one or more new supervariables were found.
+ IF (ID.GT.IDL) GO TO 210
+C Check whether every variable of element KE is in new element ME.
+ DO 200 JP = JP1,JP2
+ JS = IW(JP)
+ IF (FLAG(JS).NE.0) GO TO 210
+ 200 CONTINUE
+C Set tree pointer and FLAG to indicate that element KE is absorbed
+C into new element ME.
+ IPE(KE) = -ME
+ FLAG(KE) = -1
+ GO TO 220
+C Store element KE in the revised list for supervariable IS and flag it.
+ 210 IW(NP) = KE
+ FLAG(KE) = -NFLG
+ NP = NP + 1
+ 220 CONTINUE
+ NP0 = NP
+ GO TO 250
+C Treat the rest of the list associated with supervariable IS. It
+C consists entirely of supervariables.
+ 230 KP0 = KP
+ NP0 = NP
+ DO 240 KP = KP0,KP2
+ KS = IW(KP)
+ IF (FLAG(KS).LE.NFLG) THEN
+ IF (IPE(KS).EQ.-ROOT) THEN
+ KS = ROOT
+ IW(KP) = ROOT
+ IF (FLAG(KS).LE.NFLG) GO TO 240
+ ELSE
+ GO TO 240
+ END IF
+ END IF
+C Add to degree, flag supervariable KS and add it to new list.
+ ID = ID + NV(KS)
+ FLAG(KS) = NFLG
+ IW(NP) = KS
+ NP = NP + 1
+ 240 CONTINUE
+C Move first supervariable to end of list, move first element to end
+C of element part of list and add new element to front of list.
+ 250 IF (ID.GE.LIMIT) GO TO 295
+ IW(NP) = IW(NP0)
+ IW(NP0) = IW(KP1)
+ IW(KP1) = ME
+C Store the new length of the list.
+ IW(KP1-1) = NP - KP1 + 1
+C
+C Check whether row is is identical to another by looking in linked
+C list of supervariables with degree ID at those whose lists have
+C first entry ME. Note that those containing ME come first so the
+C search can be terminated when a list not starting with ME is
+C found.
+ JS = IPD(ID)
+ DO 280 L = 1,N
+ IF (JS.LE.0) GO TO 300
+ KP1 = IPE(JS) + 1
+ IF (IW(KP1).NE.ME) GO TO 300
+C JS has same degree and is active. Check if identical to IS.
+ KP2 = KP1 - 1 + IW(KP1-1)
+ DO 260 KP = KP1,KP2
+ IE = IW(KP)
+C Jump if IE is a supervariable or an element not in the list of IS.
+ IF (ABS(FLAG(IE)+0).GT.NFLG) GO TO 270
+ 260 CONTINUE
+ GO TO 290
+
+ 270 JS = NXT(JS)
+ 280 CONTINUE
+C Supervariable amalgamation. Row IS is identical to row JS.
+C Regard all variables in the two supervariables as being in IS. Set
+C tree pointer, FLAG and NV entries.
+ 290 IPE(JS) = -IS
+ NV(IS) = NV(IS) + NV(JS)
+ NV(JS) = 0
+ FLAG(JS) = -1
+C Replace JS by IS in linked list.
+ NS = NXT(JS)
+ LS = LST(JS)
+ IF (NS.GT.0) LST(NS) = IS
+ IF (LS.GT.0) NXT(LS) = IS
+ LST(IS) = LS
+ NXT(IS) = NS
+ LST(JS) = 0
+ NXT(JS) = 0
+ IF (IPD(ID).EQ.JS) IPD(ID) = IS
+ GO TO 310
+C Treat IS as full. Merge it into the root node.
+ 295 IF (NVROOT.EQ.0) THEN
+ ROOT = IS
+ IPE(IS) = 0
+ ELSE
+ IW(K) = ROOT
+ IPE(IS) = -ROOT
+ NV(ROOT) = NV(ROOT) + NV(IS)
+ NV(IS) = 0
+ FLAG(IS) = -1
+ END IF
+ NVROOT = NV(ROOT)
+ GO TO 310
+C Insert IS into linked list of supervariables of same degree.
+ 300 NS = IPD(ID)
+ IF (NS.GT.0) LST(NS) = IS
+ NXT(IS) = NS
+ IPD(ID) = IS
+ LST(IS) = -ID
+ MD = MIN(MD,ID)
+ 310 CONTINUE
+
+C
+C Reset flags for supervariables in newly created element and
+C remove those absorbed into others.
+ DO 320 K = K1,K2
+ IS = IW(K)
+ IF (NV(IS).EQ.0) GO TO 320
+ FLAG(IS) = NFLG
+ IW(IP) = IS
+ IP = IP + 1
+ 320 CONTINUE
+
+ FLAG(ME) = -NFLG
+C Move first entry to end to make room for length.
+ IW(IP) = IW(K1)
+ IW(K1) = IP - K1
+C Set pointer for new element and reset IWFR.
+ IPE(ME) = K1
+ IWFR = IP + 1
+
+C End of main loop
+ 340 CONTINUE
+C
+
+C Absorb any remaining variables into the root
+ 350 DO 360 IS = 1,N
+ IF(NXT(IS).NE.0 .OR. LST(IS).NE.0) THEN
+ IF (NVROOT.EQ.0) THEN
+ ROOT = IS
+ IPE(IS) = 0
+ ELSE
+ IPE(IS) = -ROOT
+ END IF
+ NVROOT = NVROOT + NV(IS)
+ NV(IS) = 0
+ END IF
+ 360 CONTINUE
+C Link any remaining elements to the root
+ DO 370 IE = 1,N
+ IF (IPE(IE).GT.0) IPE(IE) = -ROOT
+ 370 CONTINUE
+ IF(NVROOT.GT.0)NV(ROOT)=NVROOT
+ END
+ SUBROUTINE MA27UD(N,IPE,IW,LW,IWFR,NCMPA)
+C COMPRESS LISTS HELD BY MA27H/HD AND MA27K/KD IN IW AND ADJUST POINTERS
+C IN IPE TO CORRESPOND.
+C N IS THE MATRIX ORDER. IT IS NOT ALTERED.
+C IPE(I) POINTS TO THE POSITION IN IW OF THE START OF LIST I OR IS
+C ZERO IF THERE IS NO LIST I. ON EXIT IT POINTS TO THE NEW POSITION.
+C IW HOLDS THE LISTS, EACH HEADED BY ITS LENGTH. ON OUTPUT THE SAME
+C LISTS ARE HELD, BUT THEY ARE NOW COMPRESSED TOGETHER.
+C LW HOLDS THE LENGTH OF IW. IT IS NOT ALTERED.
+C IWFR NEED NOT BE SET ON ENTRY. ON EXIT IT POINTS TO THE FIRST FREE
+C LOCATION IN IW.
+C ON RETURN IT IS SET TO THE FIRST FREE LOCATION IN IW.
+C NCMPA see INFO(11) in MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER IWFR,LW,N,NCMPA
+C ..
+C .. Array Arguments ..
+ INTEGER IPE(N),IW(LW)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IR,K,K1,K2,LWFR
+C ..
+C .. Executable Statements ..
+ NCMPA = NCMPA + 1
+C PREPARE FOR COMPRESSING BY STORING THE LENGTHS OF THE
+C LISTS IN IPE AND SETTING THE FIRST ENTRY OF EACH LIST TO
+C -(LIST NUMBER).
+ DO 10 I = 1,N
+ K1 = IPE(I)
+ IF (K1.LE.0) GO TO 10
+ IPE(I) = IW(K1)
+ IW(K1) = -I
+ 10 CONTINUE
+C
+C COMPRESS
+C IWFR POINTS JUST BEYOND THE END OF THE COMPRESSED FILE.
+C LWFR POINTS JUST BEYOND THE END OF THE UNCOMPRESSED FILE.
+ IWFR = 1
+ LWFR = IWFR
+ DO 60 IR = 1,N
+ IF (LWFR.GT.LW) GO TO 70
+C SEARCH FOR THE NEXT NEGATIVE ENTRY.
+ DO 20 K = LWFR,LW
+ IF (IW(K).LT.0) GO TO 30
+ 20 CONTINUE
+ GO TO 70
+C PICK UP ENTRY NUMBER, STORE LENGTH IN NEW POSITION, SET NEW POINTER
+C AND PREPARE TO COPY LIST.
+ 30 I = -IW(K)
+ IW(IWFR) = IPE(I)
+ IPE(I) = IWFR
+ K1 = K + 1
+ K2 = K + IW(IWFR)
+ IWFR = IWFR + 1
+ IF (K1.GT.K2) GO TO 50
+C COPY LIST TO NEW POSITION.
+ DO 40 K = K1,K2
+ IW(IWFR) = IW(K)
+ IWFR = IWFR + 1
+ 40 CONTINUE
+ 50 LWFR = K2 + 1
+ 60 CONTINUE
+ 70 RETURN
+
+ END
+ SUBROUTINE MA27JD(N,NZ,IRN,ICN,PERM,IW,LW,IPE,IQ,FLAG,IWFR,
+ + ICNTL,INFO)
+C
+C SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27K/KD.
+C
+C GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC
+C MATRIX AND A PERMUTATION, CONSTRUCT THE SPARSITY PATTERN
+C OF THE STRICTLY UPPER TRIANGULAR PART OF THE PERMUTED MATRIX.
+C EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT
+C THE PAIR. DIAGONAL ELEMENTS ARE IGNORED. NO CHECK IS MADE
+C FOR DUPLICATE ELEMENTS UNLESS ANY ROW HAS MORE THAN ICNTL(4)
+C NON-ZEROS, IN WHICH CASE DUPLICATES ARE REMOVED.
+C
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT
+C ALTERED.
+C IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW INDICES OF THE
+C NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED WITH IW.
+C IRN(1) MAY BE EQUIVALENCED WITH IW(1).
+C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN INDICES OF THE
+C NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED
+C WITH IW.ICN(1) MAY BE EQUIVELENCED WITH IW(K),K.GT.NZ.
+C PERM(I) MUST BE SET TO HOLD THE POSITION OF VARIABLE I IN THE
+C PERMUTED ORDER. IT IS NOT ALTERED.
+C IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF
+C COLUMN NUMBERS, EACH LIST BEING HEADED BY ITS LENGTH.
+C LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST
+C MAX(NZ,N+(NO. OF OFF-DIAGONAL NON-ZEROS)). IT IS NOT ALTERED.
+C IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF
+C THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY.
+C IQ NEED NOT BE SET. ON OUTPUT IQ(I) CONTAINS THE NUMBER OF
+C OFF-DIAGONAL NON-ZEROS IN ROW I, INCLUDING DUPLICATES.
+C FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE
+C ENTRIES TO BE IDENTIFIED QUICKLY.
+C IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST
+C UNUSED LOCATION IN IW.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER IWFR,LW,N,NZ
+C ..
+C .. Array Arguments ..
+ INTEGER FLAG(N),ICN(*),IPE(N),IQ(N),IRN(*),IW(LW),PERM(N)
+ INTEGER ICNTL(30),INFO(20)
+C ..
+C .. Local Scalars ..
+ INTEGER I,ID,IN,J,JDUMMY,K,K1,K2,L,LBIG,LEN
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MAX
+C ..
+C .. Executable Statements ..
+C
+C INITIALIZE INFO(1), INFO(2) AND IQ
+ INFO(1) = 0
+ INFO(2) = 0
+ DO 10 I = 1,N
+ IQ(I) = 0
+ 10 CONTINUE
+C
+C COUNT THE NUMBERS OF NON-ZEROS IN THE ROWS, PRINT WARNINGS ABOUT
+C OUT-OF-RANGE INDICES AND TRANSFER GENUINE ROW NUMBERS
+C (NEGATED) INTO IW.
+ IF (NZ.EQ.0) GO TO 110
+ DO 100 K = 1,NZ
+ I = IRN(K)
+ J = ICN(K)
+ IW(K) = -I
+ IF(I.LT.J) THEN
+ IF (I.GE.1 .AND. J.LE.N) GO TO 80
+ ELSE IF(I.GT.J) THEN
+ IF (J.GE.1 .AND. I.LE.N) GO TO 80
+ ELSE
+ IW(K) = 0
+ IF (I.GE.1 .AND. I.LE.N) GO TO 100
+ END IF
+ INFO(2) = INFO(2) + 1
+ INFO(1) = 1
+ IW(K) = 0
+ IF (INFO(2).LE.1 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=60) INFO(1)
+ END IF
+
+ 60 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27AD',
+ + ' *** INFO(1) =',I2)
+
+ IF (INFO(2).LE.10 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=70) K,I,J
+ END IF
+
+ 70 FORMAT (I6,'TH NON-ZERO (IN ROW',I6,' AND COLUMN',I6,
+ + ') IGNORED')
+
+ GO TO 100
+
+ 80 IF (PERM(J).GT.PERM(I)) GO TO 90
+ IQ(J) = IQ(J) + 1
+ GO TO 100
+
+ 90 IQ(I) = IQ(I) + 1
+ 100 CONTINUE
+C
+C ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW ENDS
+C IN IPE.
+ 110 IWFR = 1
+ LBIG = 0
+ DO 120 I = 1,N
+ L = IQ(I)
+ LBIG = MAX(L,LBIG)
+ IWFR = IWFR + L
+ IPE(I) = IWFR - 1
+ 120 CONTINUE
+C
+C PERFORM IN-PLACE SORT
+ IF (NZ.EQ.0) GO TO 250
+ DO 160 K = 1,NZ
+ I = -IW(K)
+ IF (I.LE.0) GO TO 160
+ L = K
+ IW(K) = 0
+ DO 150 ID = 1,NZ
+ J = ICN(L)
+ IF (PERM(I).LT.PERM(J)) GO TO 130
+ L = IPE(J)
+ IPE(J) = L - 1
+ IN = IW(L)
+ IW(L) = I
+ GO TO 140
+
+ 130 L = IPE(I)
+ IPE(I) = L - 1
+ IN = IW(L)
+ IW(L) = J
+ 140 I = -IN
+ IF (I.LE.0) GO TO 160
+ 150 CONTINUE
+ 160 CONTINUE
+C
+C MAKE ROOM IN IW FOR ROW LENGTHS AND INITIALIZE FLAG.
+ K = IWFR - 1
+ L = K + N
+ IWFR = L + 1
+ DO 190 I = 1,N
+ FLAG(I) = 0
+ J = N + 1 - I
+ LEN = IQ(J)
+ IF (LEN.LE.0) GO TO 180
+ DO 170 JDUMMY = 1,LEN
+ IW(L) = IW(K)
+ K = K - 1
+ L = L - 1
+ 170 CONTINUE
+ 180 IPE(J) = L
+ L = L - 1
+ 190 CONTINUE
+ IF (LBIG.GE.ICNTL(4)) GO TO 210
+C
+C PLACE ROW LENGTHS IN IW
+ DO 200 I = 1,N
+ K = IPE(I)
+ IW(K) = IQ(I)
+ IF (IQ(I).EQ.0) IPE(I) = 0
+ 200 CONTINUE
+ GO TO 250
+C
+C
+C REMOVE DUPLICATE ENTRIES
+ 210 IWFR = 1
+ DO 240 I = 1,N
+ K1 = IPE(I) + 1
+ K2 = IPE(I) + IQ(I)
+ IF (K1.LE.K2) GO TO 220
+ IPE(I) = 0
+ GO TO 240
+
+ 220 IPE(I) = IWFR
+ IWFR = IWFR + 1
+ DO 230 K = K1,K2
+ J = IW(K)
+ IF (FLAG(J).EQ.I) GO TO 230
+ IW(IWFR) = J
+ IWFR = IWFR + 1
+ FLAG(J) = I
+ 230 CONTINUE
+ K = IPE(I)
+ IW(K) = IWFR - K - 1
+ 240 CONTINUE
+ 250 RETURN
+
+ END
+ SUBROUTINE MA27KD(N,IPE,IW,LW,IWFR,IPS,IPV,NV,FLAG,NCMPA)
+C
+C USING A GIVEN PIVOTAL SEQUENCE AND A REPRESENTATION OF THE MATRIX THAT
+C INCLUDES ONLY NON-ZEROS OF THE STRICTLY UPPER-TRIANGULAR PART
+C OF THE PERMUTED MATRIX, CONSTRUCT TREE POINTERS.
+C
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE
+C START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF-
+C DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS.
+C IF VARIABLE I IS ELIMINATED THEN IPE(I) POINTS TO THE LIST
+C OF VARIABLES FOR CREATED ELEMENT I. IF ELEMENT I IS
+C ABSORBED INTO NEWLY CREATED ELEMENT J THEN IPE(I)=-J.
+C IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY
+C ROWS, EACH LIST BEING HEADED BY ITS LENGTH. WHEN A VARIABLE
+C IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF VARIABLES
+C IN THE NEW ELEMENT.
+C LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED.
+C IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE.
+C IT IS REVISED DURING EXECUTION, CONTINUING TO HAVE THIS MEANING.
+C IPS(I) MUST BE SET TO THE POSITION OF VARIABLE I IN THE REQUIRED
+C ORDERING. IT IS NOT ALTERED.
+C IPV NEED NOT BE SET. IPV(K) IS SET TO HOLD THE K TH VARIABLE
+C IN PIVOT ORDER.
+C NV NEED NOT BE SET. IF VARIABLE J HAS NOT BEEN ELIMINATED THEN
+C THE LAST ELEMENT WHOSE LEADING VARIABLE (VARIABLE EARLIEST
+C IN THE PIVOT SEQUENCE) IS J IS ELEMENT NV(J). IF ELEMENT J
+C EXISTS THEN THE LAST ELEMENT HAVING THE SAME LEADING
+C VARIABLE IS NV(J). IN BOTH CASES NV(J)=0 IF THERE IS NO SUCH
+C ELEMENT. IF ELEMENT J HAS BEEN MERGED INTO A LATER ELEMENT
+C THEN NV(J) IS THE DEGREE AT THE TIME OF ELIMINATION.
+C FLAG IS USED AS WORKSPACE FOR VARIABLE FLAGS.
+C FLAG(JS)=ME IF JS HAS BEEN INCLUDED IN THE LIST FOR ME.
+C NCMPA see INFO(11) in MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER IWFR,LW,N,NCMPA
+C ..
+C .. Array Arguments ..
+ INTEGER FLAG(N),IPE(N),IPS(N),IPV(N),IW(LW),NV(N)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IE,IP,J,JE,JP,JP1,JP2,JS,KDUMMY,LN,LWFR,ME,MINJS,ML,MS
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27UD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MIN
+C ..
+C .. Executable Statements ..
+C
+C INITIALIZATIONS
+ DO 10 I = 1,N
+ FLAG(I) = 0
+ NV(I) = 0
+ J = IPS(I)
+ IPV(J) = I
+ 10 CONTINUE
+ NCMPA = 0
+C
+C START OF MAIN LOOP
+C
+ DO 100 ML = 1,N
+C ME=MS IS THE NAME OF THE VARIABLE ELIMINATED AND
+C OF THE ELEMENT CREATED IN THE MAIN LOOP.
+ MS = IPV(ML)
+ ME = MS
+ FLAG(MS) = ME
+C
+C MERGE ROW MS WITH ALL THE ELEMENTS HAVING MS AS LEADING VARIABLE.
+C IP POINTS TO THE START OF THE NEW LIST.
+ IP = IWFR
+C MINJS IS SET TO THE POSITION IN THE ORDER OF THE LEADING VARIABLE
+C IN THE NEW LIST.
+ MINJS = N
+ IE = ME
+ DO 70 KDUMMY = 1,N
+C SEARCH VARIABLE LIST OF ELEMENT IE.
+C JP POINTS TO THE CURRENT POSITION IN THE LIST BEING SEARCHED.
+ JP = IPE(IE)
+C LN IS THE LENGTH OF THE LIST BEING SEARCHED.
+ LN = 0
+ IF (JP.LE.0) GO TO 60
+ LN = IW(JP)
+C
+C SEARCH FOR DIFFERENT VARIABLES AND ADD THEM TO LIST,
+C COMPRESSING WHEN NECESSARY
+ DO 50 JP1 = 1,LN
+ JP = JP + 1
+C PLACE NEXT VARIABLE IN JS.
+ JS = IW(JP)
+C JUMP IF VARIABLE HAS ALREADY BEEN INCLUDED.
+ IF (FLAG(JS).EQ.ME) GO TO 50
+ FLAG(JS) = ME
+ IF (IWFR.LT.LW) GO TO 40
+C PREPARE FOR COMPRESSING IW BY ADJUSTING POINTER TO AND LENGTH OF
+C THE LIST FOR IE TO REFER TO THE REMAINING ENTRIES.
+ IPE(IE) = JP
+ IW(JP) = LN - JP1
+C COMPRESS IW.
+ CALL MA27UD(N,IPE,IW,IP-1,LWFR,NCMPA)
+C COPY NEW LIST FORWARD
+ JP2 = IWFR - 1
+ IWFR = LWFR
+ IF (IP.GT.JP2) GO TO 30
+ DO 20 JP = IP,JP2
+ IW(IWFR) = IW(JP)
+ IWFR = IWFR + 1
+ 20 CONTINUE
+ 30 IP = LWFR
+ JP = IPE(IE)
+C ADD VARIABLE JS TO NEW LIST.
+ 40 IW(IWFR) = JS
+ MINJS = MIN(MINJS,IPS(JS)+0)
+ IWFR = IWFR + 1
+ 50 CONTINUE
+C RECORD ABSORPTION OF ELEMENT IE INTO NEW ELEMENT.
+ 60 IPE(IE) = -ME
+C PICK UP NEXT ELEMENT WITH LEADING VARIABLE MS.
+ JE = NV(IE)
+C STORE DEGREE OF IE.
+ NV(IE) = LN + 1
+ IE = JE
+C LEAVE LOOP IF THERE ARE NO MORE ELEMENTS.
+ IF (IE.EQ.0) GO TO 80
+ 70 CONTINUE
+ 80 IF (IWFR.GT.IP) GO TO 90
+C DEAL WITH NULL NEW ELEMENT.
+ IPE(ME) = 0
+ NV(ME) = 1
+ GO TO 100
+C LINK NEW ELEMENT WITH OTHERS HAVING SAME LEADING VARIABLE.
+ 90 MINJS = IPV(MINJS)
+ NV(ME) = NV(MINJS)
+ NV(MINJS) = ME
+C MOVE FIRST ENTRY IN NEW LIST TO END TO ALLOW ROOM FOR LENGTH AT
+C FRONT. SET POINTER TO FRONT.
+ IW(IWFR) = IW(IP)
+ IW(IP) = IWFR - IP
+ IPE(ME) = IP
+ IWFR = IWFR + 1
+ 100 CONTINUE
+ RETURN
+
+ END
+ SUBROUTINE MA27LD(N,IPE,NV,IPS,NE,NA,ND,NSTEPS,NEMIN)
+C
+C TREE SEARCH
+C
+C GIVEN SON TO FATHER TREE POINTERS, PERFORM DEPTH-FIRST
+C SEARCH TO FIND PIVOT ORDER AND NUMBER OF ELIMINATIONS
+C AND ASSEMBLIES AT EACH STAGE.
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C IPE(I) MUST BE SET EQUAL TO -(FATHER OF NODE I) OR ZERO IF
+C NODE IS A ROOT. IT IS ALTERED TO POINT TO ITS NEXT
+C YOUNGER BROTHER IF IT HAS ONE, BUT OTHERWISE IS NOT
+C CHANGED.
+C NV(I) MUST BE SET TO ZERO IF NO VARIABLES ARE ELIMINATED AT NODE
+C I AND TO THE DEGREE OTHERWISE. ONLY LEAF NODES CAN HAVE
+C ZERO VALUES OF NV(I). NV IS NOT ALTERED.
+C IPS(I) NEED NOT BE SET. IT IS USED TEMPORARILY TO HOLD
+C -(ELDEST SON OF NODE I) IF IT HAS ONE AND 0 OTHERWISE. IT IS
+C EVENTUALLY SET TO HOLD THE POSITION OF NODE I IN THE ORDER.
+C NE(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF VARIABLES
+C ELIMINATED AT STAGE IS OF THE ELIMINATION.
+C NA(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELEMENTS
+C ASSEMBLED AT STAGE IS OF THE ELIMINATION.
+C ND(IS) NEED NOT BE SET. IT IS SET TO THE DEGREE AT STAGE IS OF
+C THE ELIMINATION.
+C NSTEPS NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELIMINATION
+C STEPS.
+C NEMIN see ICNTL(5) in MA27A/AD.
+C
+C .. Scalar Arguments ..
+ INTEGER N,NSTEPS,NEMIN
+C ..
+C .. Array Arguments ..
+ INTEGER IPE(N),IPS(N),NA(N),ND(N),NE(N),NV(N)
+C ..
+C .. Local Scalars ..
+ INTEGER I,IB,IF,IL,IS,ISON,K,L,NR
+C ..
+C .. Executable Statements ..
+C INITIALIZE IPS AND NE.
+ DO 10 I = 1,N
+ IPS(I) = 0
+ NE(I) = 0
+ 10 CONTINUE
+C
+C SET IPS(I) TO -(ELDEST SON OF NODE I) AND IPE(I) TO NEXT YOUNGER
+C BROTHER OF NODE I IF IT HAS ONE.
+C FIRST PASS IS FOR NODES WITHOUT ELIMINATIONS.
+ DO 20 I = 1,N
+ IF (NV(I).GT.0) GO TO 20
+ IF = -IPE(I)
+ IS = -IPS(IF)
+ IF (IS.GT.0) IPE(I) = IS
+ IPS(IF) = -I
+ 20 CONTINUE
+C NR IS DECREMENTED FOR EACH ROOT NODE. THESE ARE STORED IN
+C NE(I),I=NR,N.
+ NR = N + 1
+C SECOND PASS TO ADD NODES WITH ELIMINATIONS.
+ DO 50 I = 1,N
+ IF (NV(I).LE.0) GO TO 50
+C NODE IF IS THE FATHER OF NODE I.
+ IF = -IPE(I)
+ IF (IF.EQ.0) GO TO 40
+ IS = -IPS(IF)
+C JUMP IF NODE IF HAS NO SONS YET.
+ IF (IS.LE.0) GO TO 30
+C SET POINTER TO NEXT BROTHER
+ IPE(I) = IS
+C NODE I IS ELDEST SON OF NODE IF.
+ 30 IPS(IF) = -I
+ GO TO 50
+C WE HAVE A ROOT
+ 40 NR = NR - 1
+ NE(NR) = I
+ 50 CONTINUE
+C
+C DEPTH-FIRST SEARCH.
+C IL HOLDS THE CURRENT TREE LEVEL. ROOTS ARE AT LEVEL N, THEIR SONS
+C ARE AT LEVEL N-1, ETC.
+C IS HOLDS THE CURRENT ELIMINATION STAGE. WE ACCUMULATE THE NUMBER
+C OF ELIMINATIONS AT STAGE IS DIRECTLY IN NE(IS). THE NUMBER OF
+C ASSEMBLIES IS ACCUMULATED TEMPORARILY IN NA(IL), FOR TREE
+C LEVEL IL, AND IS TRANSFERED TO NA(IS) WHEN WE REACH THE
+C APPROPRIATE STAGE IS.
+ IS = 1
+C I IS THE CURRENT NODE.
+ I = 0
+ DO 160 K = 1,N
+ IF (I.GT.0) GO TO 60
+C PICK UP NEXT ROOT.
+ I = NE(NR)
+ NE(NR) = 0
+ NR = NR + 1
+ IL = N
+ NA(N) = 0
+C GO TO SON FOR AS LONG AS POSSIBLE, CLEARING FATHER-SON POINTERS
+C IN IPS AS EACH IS USED AND SETTING NA(IL)=0 FOR ALL LEVELS
+C REACHED.
+ 60 DO 70 L = 1,N
+ IF (IPS(I).GE.0) GO TO 80
+ ISON = -IPS(I)
+ IPS(I) = 0
+ I = ISON
+ IL = IL - 1
+ NA(IL) = 0
+ 70 CONTINUE
+C RECORD POSITION OF NODE I IN THE ORDER.
+ 80 IPS(I) = K
+ NE(IS) = NE(IS) + 1
+C JUMP IF NODE HAS NO ELIMINATIONS.
+ IF (NV(I).LE.0) GO TO 120
+ IF (IL.LT.N) NA(IL+1) = NA(IL+1) + 1
+ NA(IS) = NA(IL)
+ ND(IS) = NV(I)
+C CHECK FOR STATIC CONDENSATION
+ IF (NA(IS).NE.1) GO TO 90
+ IF (ND(IS-1)-NE(IS-1).EQ.ND(IS)) GO TO 100
+C CHECK FOR SMALL NUMBERS OF ELIMINATIONS IN BOTH LAST TWO STEPS.
+ 90 IF (NE(IS).GE.NEMIN) GO TO 110
+ IF (NA(IS).EQ.0) GO TO 110
+ IF (NE(IS-1).GE.NEMIN) GO TO 110
+C COMBINE THE LAST TWO STEPS
+ 100 NA(IS-1) = NA(IS-1) + NA(IS) - 1
+ ND(IS-1) = ND(IS) + NE(IS-1)
+ NE(IS-1) = NE(IS) + NE(IS-1)
+ NE(IS) = 0
+ GO TO 120
+
+ 110 IS = IS + 1
+ 120 IB = IPE(I)
+ IF (IB.GE.0) THEN
+C NODE I HAS A BROTHER OR IS A ROOT
+ IF (IB.GT.0) NA(IL) = 0
+ I = IB
+ ELSE
+C GO TO FATHER OF NODE I
+ I = -IB
+ IL = IL + 1
+ END IF
+ 160 CONTINUE
+ NSTEPS = IS - 1
+ RETURN
+
+ END
+ SUBROUTINE MA27MD(N,NZ,IRN,ICN,PERM,NA,NE,ND,NSTEPS,LSTKI,LSTKR,
+ + IW,INFO,OPS)
+C
+C STORAGE AND OPERATION COUNT EVALUATION.
+C
+C EVALUATE NUMBER OF OPERATIONS AND SPACE REQUIRED BY FACTORIZATION
+C USING MA27B/BD. THE VALUES GIVEN ARE EXACT ONLY IF NO NUMERICAL
+C PIVOTING IS PERFORMED AND THEN ONLY IF IRN(1) WAS NOT
+C EQUIVALENCED TO IW(1) BY THE USER BEFORE CALLING MA27A/AD. IF
+C THE EQUIVALENCE HAS BEEN MADE ONLY AN UPPER BOUND FOR NIRNEC
+C AND NRLNEC CAN BE CALCULATED ALTHOUGH THE OTHER COUNTS WILL
+C STILL BE EXACT.
+C
+C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
+C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT ALTERED.
+C IRN,ICN. UNLESS IRN(1) HAS BEEN EQUIVALENCED TO IW(1)
+C IRN,ICN MUST BE SET TO THE ROW AND COLUMN INDICES OF THE
+C NON-ZEROS ON INPUT. THEY ARE NOT ALTERED BY MA27M/MD.
+C PERM MUST BE SET TO THE POSITION IN THE PIVOT ORDER OF EACH ROW.
+C IT IS NOT ALTERED.
+C NA,NE,ND MUST BE SET TO HOLD, FOR EACH TREE NODE, THE NUMBER OF STACK
+C ELEMENTS ASSEMBLED, THE NUMBER OF ELIMINATIONS AND THE SIZE OF
+C THE ASSEMBLED FRONT MATRIX RESPECTIVELY. THEY ARE NOT ALTERED.
+C NSTEPS MUST BE SET TO HOLD THE NUMBER OF TREE NODES. IT IS NOT
+C ALTERED.
+C LSTKI IS USED AS A WORK ARRAY BY MA27M/MD.
+C LSTKR. IF IRN(1) IS EQUIVALENCED TO IW(1) THEN LSTKR(I)
+C MUST HOLD THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES (INCLUDING
+C DUPLICATES) IN ROW I (I=1,..,N) OF THE ORIGINAL MATRIX. IT
+C IS USED AS WORKSPACE BY MA27M/MD.
+C IW IS A WORKSPACE ARRAY USED BY OTHER SUBROUTINES AND PASSED TO THIS
+C SUBROUTINE ONLY SO THAT A TEST FOR EQUIVALENCE WITH IRN CAN BE
+C MADE.
+C
+C COUNTS FOR OPERATIONS AND STORAGE ARE ACCUMULATED IN VARIABLES
+C OPS,NRLTOT,NIRTOT,NRLNEC,NIRNEC,NRLADU,NRLNEC,NIRNEC.
+C OPS NUMBER OF MULTIPLICATIONS AND ADDITIONS DURING FACTORIZATION.
+C NRLADU,NIRADU REAL AND INTEGER STORAGE RESPECTIVELY FOR THE
+C MATRIX FACTORS.
+C NRLTOT,NIRTOT REAL AND INTEGER STRORAGE RESPECTIVELY REQUIRED
+C FOR THE FACTORIZATION IF NO COMPRESSES ARE ALLOWED.
+C NRLNEC,NIRNEC REAL AND INTEGER STORAGE RESPECTIVELY REQUIRED FOR
+C THE FACTORIZATION IF COMPRESSES ARE ALLOWED.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE
+C TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE.
+C
+C .. Scalar Arguments ..
+ DOUBLE PRECISION OPS
+ INTEGER N,NSTEPS,NZ
+C ..
+C .. Array Arguments ..
+ INTEGER ICN(*),IRN(*),IW(*),LSTKI(N),LSTKR(N),NA(NSTEPS),
+ + ND(NSTEPS),NE(NSTEPS),PERM(N),INFO(20)
+C ..
+C .. Local Scalars ..
+ INTEGER I,INEW,IOLD,IORG,IROW,ISTKI,ISTKR,ITOP,ITREE,JOLD,JORG,K,
+ + LSTK,NASSR,NELIM,NFR,NSTK,NUMORG,NZ1,NZ2
+ DOUBLE PRECISION DELIM
+ INTEGER NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MAX,MIN
+C ..
+C .. Executable Statements ..
+C
+ IF (NZ.EQ.0) GO TO 20
+C JUMP IF IW AND IRN HAVE NOT BEEN EQUIVALENCED.
+ IF (IRN(1).NE.IW(1)) GO TO 20
+C RESET IRN(1).
+ IRN(1) = IW(1) - 1
+C THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES IS ACCUMULATED IN NZ2.
+C LSTKI IS SET TO HOLD THE TOTAL NUMBER OF ENTRIES (INCUDING
+C THE DIAGONAL) IN EACH ROW IN PERMUTED ORDER.
+ NZ2 = 0
+ DO 10 IOLD = 1,N
+ INEW = PERM(IOLD)
+ LSTKI(INEW) = LSTKR(IOLD) + 1
+ NZ2 = NZ2 + LSTKR(IOLD)
+ 10 CONTINUE
+C NZ1 IS THE NUMBER OF ENTRIES IN ONE TRIANGLE INCLUDING THE DIAGONAL.
+C NZ2 IS THE TOTAL NUMBER OF ENTRIES INCLUDING THE DIAGONAL.
+ NZ1 = NZ2/2 + N
+ NZ2 = NZ2 + N
+ GO TO 60
+C COUNT (IN LSTKI) NON-ZEROS IN ORIGINAL MATRIX BY PERMUTED ROW (UPPER
+C TRIANGLE ONLY). INITIALIZE COUNTS.
+ 20 DO 30 I = 1,N
+ LSTKI(I) = 1
+ 30 CONTINUE
+C ACCUMULATE NUMBER OF NON-ZEROS WITH INDICES IN RANGE IN NZ1
+C DUPLICATES ON THE DIAGONAL ARE IGNORED BUT NZ1 INCLUDES ANY
+C DIAGONALS NOT PRESENT ON INPUT.
+C ACCUMULATE ROW COUNTS IN LSTKI.
+ NZ1 = N
+ IF (NZ.EQ.0) GO TO 50
+ DO 40 I = 1,NZ
+ IOLD = IRN(I)
+ JOLD = ICN(I)
+C JUMP IF INDEX IS OUT OF RANGE.
+ IF (IOLD.LT.1 .OR. IOLD.GT.N) GO TO 40
+ IF (JOLD.LT.1 .OR. JOLD.GT.N) GO TO 40
+ IF (IOLD.EQ.JOLD) GO TO 40
+ NZ1 = NZ1 + 1
+ IROW = MIN(PERM(IOLD)+0,PERM(JOLD)+0)
+ LSTKI(IROW) = LSTKI(IROW) + 1
+ 40 CONTINUE
+ 50 NZ2 = NZ1
+C ISTKR,ISTKI CURRENT NUMBER OF STACK ENTRIES IN
+C REAL AND INTEGER STORAGE RESPECTIVELY.
+C OPS,NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC,NZ2 ARE DEFINED ABOVE.
+C NZ2 CURRENT NUMBER OF ORIGINAL MATRIX ENTRIES NOT YET PROCESSED.
+C NUMORG CURRENT TOTAL NUMBER OF ROWS ELIMINATED.
+C ITOP CURRENT NUMBER OF ELEMENTS ON THE STACK.
+ 60 ISTKI = 0
+ ISTKR = 0
+ OPS = 0.0D0
+ NRLADU = 0
+C ONE LOCATION IS NEEDED TO RECORD THE NUMBER OF BLOCKS
+C ACTUALLY USED.
+ NIRADU = 1
+ NIRTOT = NZ1
+ NRLTOT = NZ1
+ NIRNEC = NZ2
+ NRLNEC = NZ2
+ NUMORG = 0
+ ITOP = 0
+C
+C EACH PASS THROUGH THIS LOOP PROCESSES A NODE OF THE TREE.
+ DO 100 ITREE = 1,NSTEPS
+ NELIM = NE(ITREE)
+ DELIM = NELIM
+ NFR = ND(ITREE)
+ NSTK = NA(ITREE)
+C ADJUST STORAGE COUNTS ON ASSEMBLY OF CURRENT FRONTAL MATRIX.
+ NASSR = NFR* (NFR+1)/2
+ IF (NSTK.NE.0) NASSR = NASSR - LSTKR(ITOP) + 1
+ NRLTOT = MAX(NRLTOT,NRLADU+NASSR+ISTKR+NZ1)
+ NIRTOT = MAX(NIRTOT,NIRADU+NFR+2+ISTKI+NZ1)
+ NRLNEC = MAX(NRLNEC,NRLADU+NASSR+ISTKR+NZ2)
+ NIRNEC = MAX(NIRNEC,NIRADU+NFR+2+ISTKI+NZ2)
+C DECREASE NZ2 BY THE NUMBER OF ENTRIES IN ROWS BEING ELIMINATED AT
+C THIS STAGE.
+ DO 70 IORG = 1,NELIM
+ JORG = NUMORG + IORG
+ NZ2 = NZ2 - LSTKI(JORG)
+ 70 CONTINUE
+ NUMORG = NUMORG + NELIM
+C JUMP IF THERE ARE NO STACK ASSEMBLIES AT THIS NODE.
+ IF (NSTK.LE.0) GO TO 90
+C REMOVE ELEMENTS FROM THE STACK. THERE ARE ITOP ELEMENTS ON THE
+C STACK WITH THE APPROPRIATE ENTRIES IN LSTKR,LSTKI GIVING
+C THE REAL AND INTEGER STORAGE RESPECTIVELY FOR EACH STACK
+C ELEMENT.
+ DO 80 K = 1,NSTK
+ LSTK = LSTKR(ITOP)
+ ISTKR = ISTKR - LSTK
+ LSTK = LSTKI(ITOP)
+ ISTKI = ISTKI - LSTK
+ ITOP = ITOP - 1
+ 80 CONTINUE
+C ACCUMULATE NON-ZEROS IN FACTORS AND NUMBER OF OPERATIONS.
+ 90 NRLADU = NRLADU + (NELIM* (2*NFR-NELIM+1))/2
+ NIRADU = NIRADU + 2 + NFR
+ IF (NELIM.EQ.1) NIRADU = NIRADU - 1
+ OPS = OPS + ((NFR*DELIM*(NFR+1)-(2*NFR+1)*DELIM*(DELIM+1)/2+
+ + DELIM* (DELIM+1)* (2*DELIM+1)/6)/2)
+ IF (ITREE.EQ.NSTEPS) GO TO 100
+C JUMP IF ALL OF FRONTAL MATRIX HAS BEEN ELIMINATED.
+ IF (NFR.EQ.NELIM) GO TO 100
+C STACK REMAINDER OF ELEMENT.
+ ITOP = ITOP + 1
+ LSTKR(ITOP) = (NFR-NELIM)* (NFR-NELIM+1)/2
+ LSTKI(ITOP) = NFR - NELIM + 1
+ ISTKI = ISTKI + LSTKI(ITOP)
+ ISTKR = ISTKR + LSTKR(ITOP)
+C WE DO NOT NEED TO ADJUST THE COUNTS FOR THE REAL STORAGE BECAUSE
+C THE REMAINDER OF THE FRONTAL MATRIX IS SIMPLY MOVED IN THE
+C STORAGE FROM FACTORS TO STACK AND NO EXTRA STORAGE IS REQUIRED.
+ NIRTOT = MAX(NIRTOT,NIRADU+ISTKI+NZ1)
+ NIRNEC = MAX(NIRNEC,NIRADU+ISTKI+NZ2)
+ 100 CONTINUE
+C
+C ADJUST THE STORAGE COUNTS TO ALLOW FOR THE USE OF THE REAL AND
+C INTEGER STORAGE FOR PURPOSES OTHER THAN PURELY THE
+C FACTORIZATION ITSELF.
+C THE SECOND TWO TERMS ARE THE MINUMUM AMOUNT REQUIRED BY MA27N/ND.
+ NRLNEC = MAX(NRLNEC,N+MAX(NZ,NZ1))
+ NRLTOT = MAX(NRLTOT,N+MAX(NZ,NZ1))
+ NRLNEC = MIN(NRLNEC,NRLTOT)
+ NIRNEC = MAX(NZ,NIRNEC)
+ NIRTOT = MAX(NZ,NIRTOT)
+ NIRNEC = MIN(NIRNEC,NIRTOT)
+
+ INFO(3) = NRLTOT
+ INFO(4) = NIRTOT
+ INFO(5) = NRLNEC
+ INFO(6) = NIRNEC
+ INFO(7) = NRLADU
+ INFO(8) = NIRADU
+ RETURN
+
+ END
+ SUBROUTINE MA27ND(N,NZ,NZ1,A,LA,IRN,ICN,IW,LIW,PERM,IW2,ICNTL,
+ + INFO)
+C
+C SORT PRIOR TO FACTORIZATION USING MA27O/OD.
+C
+C THIS SUBROUTINE REORDERS THE USER'S INPUT SO THAT THE UPPER TRIANGLE
+C OF THE PERMUTED MATRIX, INCLUDING THE DIAGONAL, IS
+C HELD ORDERED BY ROWS AT THE END OF THE STORAGE FOR A AND IW.
+C IT IGNORES ENTRIES WITH ONE OR BOTH INDICES OUT OF RANGE AND
+C ACCUMULATES DIAGONAL ENTRIES.
+C IT ADDS EXPLICIT ZEROS ON THE DIAGONAL WHERE NECESSARY.
+C N - MUST BE SET TO THE ORDER OF THE MATRIX.
+C IT IS NOT ALTERED BY MA27N/ND.
+C NZ - ON ENTRY NZ MUST BE SET TO THE NUMBER
+C OF NON-ZEROS INPUT. NOT ALTERED BY MA27N/ND.
+C NZ1 - ON EXIT NZ1 WILL BE EQUAL TO THE NUMBER OF ENTRIES IN THE
+C SORTED MATRIX.
+C A - ON ENTRY A(I) MUST
+C HOLD THE VALUE OF THE ORIGINAL MATRIX ELEMENT IN POSITION
+C (IRN(I),ICN(I)),I=1,NZ. ON EXIT A(LA-NZ1+I),I=1,NZ1 HOLDS
+C THE UPPER TRIANGLE OF THE PERMUTED MATRIX BY ROWS WITH
+C THE DIAGONAL ENTRY FIRST ALTHOUGH THERE IS NO FURTHER
+C ORDERING WITHIN THE ROWS THEMSELVES.
+C LA - LENGTH OF ARRAY A. MUST BE AT LEAST N+MAX(NZ,NZ1).
+C IT IS NOT ALTERED BY MA27N/ND.
+C IRN - IRN(I) MUST BE SET TO
+C HOLD THE ROW INDEX OF ENTRY A(I),I=1,NZ. IRN WILL BE
+C UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW.
+C ICN - ICN(I) MUST BE SET TO
+C HOLD THE COLUMN INDEX OF ENTRY A(I),I=1,NZ. ICN WILL BE
+C UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW.
+C IW - USED AS WORKSPACE AND ON
+C EXIT, ENTRIES IW(LIW-NZ1+I),I=1,NZ1 HOLD THE COLUMN INDICES
+C (THE ORIGINAL UNPERMUTED INDICES) OF THE CORRESPONDING ENTRY
+C OF A WITH THE FIRST ENTRY FOR EACH ROW FLAGGED NEGATIVE.
+C IRN(1) MAY BE EQUIVALENCED WITH IW(1) AND ICN(1) MAY BE
+C EQUIVALENCED WITH IW(K) WHERE K.GT.NZ.
+C LIW - LENGTH OF ARRAY IW. MUST BE AT LEAST AS
+C GREAT AS THE MAXIMUM OF NZ AND NZ1.
+C NOT ALTERED BY MA27N/ND.
+C PERM - PERM(I) HOLDS THE
+C POSITION IN THE TENTATIVE PIVOT ORDER OF ROW I IN THE
+C ORIGINAL MATRIX,I=1,N. NOT ALTERED BY MA27N/ND.
+C IW2 - USED AS WORKSPACE.
+C SEE COMMENTS IN CODE IMMEDIATELY PRIOR TO
+C EACH USE.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C INFO(1) - ON EXIT FROM MA27N/ND, A ZERO VALUE OF
+C INFO(1) INDICATES THAT NO ERROR HAS BEEN DETECTED.
+C POSSIBLE NON-ZERO VALUES ARE ..
+C +1 WARNING. INDICES OUT OF RANGE. THESE ARE IGNORED,
+C THEIR NUMBER IS RECORDED IN INFO(2) OF MA27E/ED AND
+C MESSAGES IDENTIFYING THE FIRST TEN ARE OUTPUT ON UNIT
+C ICNTL(2).
+C -3 INTEGER ARRAY IW IS TOO SMALL.
+C -4 DOUBLE PRECISION ARRAY A IS TOO SMALL.
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER (ZERO=0.0D0)
+C ..
+C .. Scalar Arguments ..
+ INTEGER LA,LIW,N,NZ,NZ1
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA)
+ INTEGER ICN(*),IRN(*),IW(LIW),IW2(N),PERM(N),ICNTL(30),INFO(20)
+C ..
+C .. Local Scalars ..
+ DOUBLE PRECISION ANEXT,ANOW
+ INTEGER I,IA,ICH,II,IIW,INEW,IOLD,IPOS,J1,J2,JJ,JNEW,JOLD,JPOS,K
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MIN
+C ..
+C .. Executable Statements ..
+ INFO(1) = 0
+C INITIALIZE WORK ARRAY (IW2) IN PREPARATION FOR
+C COUNTING NUMBERS OF NON-ZEROS IN THE ROWS AND INITIALIZE
+C LAST N ENTRIES IN A WHICH WILL HOLD THE DIAGONAL ENTRIES
+ IA = LA
+ DO 10 IOLD = 1,N
+ IW2(IOLD) = 1
+ A(IA) = ZERO
+ IA = IA - 1
+ 10 CONTINUE
+C SCAN INPUT COPYING ROW INDICES FROM IRN TO THE FIRST NZ POSITIONS
+C IN IW. THE NEGATIVE OF THE INDEX IS HELD TO FLAG ENTRIES FOR
+C THE IN-PLACE SORT. ENTRIES IN IW CORRESPONDING TO DIAGONALS AND
+C ENTRIES WITH OUT-OF-RANGE INDICES ARE SET TO ZERO.
+C FOR DIAGONAL ENTRIES, REALS ARE ACCUMULATED IN THE LAST N
+C LOCATIONS OF A.
+C THE NUMBER OF ENTRIES IN EACH ROW OF THE PERMUTED MATRIX IS
+C ACCUMULATED IN IW2.
+C INDICES OUT OF RANGE ARE IGNORED AFTER BEING COUNTED AND
+C AFTER APPROPRIATE MESSAGES HAVE BEEN PRINTED.
+ INFO(2) = 0
+C NZ1 IS THE NUMBER OF NON-ZEROS HELD AFTER INDICES OUT OF RANGE HAVE
+C BEEN IGNORED AND DIAGONAL ENTRIES ACCUMULATED.
+ NZ1 = N
+ IF (NZ.EQ.0) GO TO 80
+ DO 70 K = 1,NZ
+ IOLD = IRN(K)
+ IF (IOLD.GT.N .OR. IOLD.LE.0) GO TO 30
+ JOLD = ICN(K)
+ IF (JOLD.GT.N .OR. JOLD.LE.0) GO TO 30
+ INEW = PERM(IOLD)
+ JNEW = PERM(JOLD)
+ IF (INEW.NE.JNEW) GO TO 20
+ IA = LA - N + IOLD
+ A(IA) = A(IA) + A(K)
+ GO TO 60
+
+ 20 INEW = MIN(INEW,JNEW)
+C INCREMENT NUMBER OF ENTRIES IN ROW INEW.
+ IW2(INEW) = IW2(INEW) + 1
+ IW(K) = -IOLD
+ NZ1 = NZ1 + 1
+ GO TO 70
+C ENTRY OUT OF RANGE. IT WILL BE IGNORED AND A FLAG SET.
+ 30 INFO(1) = 1
+ INFO(2) = INFO(2) + 1
+ IF (INFO(2).LE.1 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=40) INFO(1)
+ ENDIF
+
+ 40 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
+ + ' *** INFO(1) =',I2)
+
+ IF (INFO(2).LE.10 .AND. ICNTL(2).GT.0) THEN
+ WRITE (ICNTL(2),FMT=50) K,IRN(K),ICN(K)
+ END IF
+
+ 50 FORMAT (I6,'TH NON-ZERO (IN ROW',I6,' AND COLUMN',I6,
+ + ') IGNORED')
+
+ 60 IW(K) = 0
+ 70 CONTINUE
+C CALCULATE POINTERS (IN IW2) TO THE POSITION IMMEDIATELY AFTER THE END
+C OF EACH ROW.
+ 80 IF (NZ.LT.NZ1 .AND. NZ1.NE.N) GO TO 100
+C ROOM IS INCLUDED FOR THE DIAGONALS.
+ K = 1
+ DO 90 I = 1,N
+ K = K + IW2(I)
+ IW2(I) = K
+ 90 CONTINUE
+ GO TO 120
+C ROOM IS NOT INCLUDED FOR THE DIAGONALS.
+ 100 K = 1
+ DO 110 I = 1,N
+ K = K + IW2(I) - 1
+ IW2(I) = K
+ 110 CONTINUE
+C FAIL IF INSUFFICIENT SPACE IN ARRAYS A OR IW.
+ 120 IF (NZ1.GT.LIW) GO TO 210
+ IF (NZ1+N.GT.LA) GO TO 220
+C NOW RUN THROUGH NON-ZEROS IN ORDER PLACING THEM IN THEIR NEW
+C POSITION AND DECREMENTING APPROPRIATE IW2 ENTRY. IF WE ARE
+C ABOUT TO OVERWRITE AN ENTRY NOT YET MOVED, WE MUST DEAL WITH
+C THIS AT THIS TIME.
+ IF (NZ1.EQ.N) GO TO 180
+ DO 140 K = 1,NZ
+ IOLD = -IW(K)
+ IF (IOLD.LE.0) GO TO 140
+ JOLD = ICN(K)
+ ANOW = A(K)
+ IW(K) = 0
+ DO 130 ICH = 1,NZ
+ INEW = PERM(IOLD)
+ JNEW = PERM(JOLD)
+ INEW = MIN(INEW,JNEW)
+ IF (INEW.EQ.PERM(JOLD)) JOLD = IOLD
+ JPOS = IW2(INEW) - 1
+ IOLD = -IW(JPOS)
+ ANEXT = A(JPOS)
+ A(JPOS) = ANOW
+ IW(JPOS) = JOLD
+ IW2(INEW) = JPOS
+ IF (IOLD.EQ.0) GO TO 140
+ ANOW = ANEXT
+ JOLD = ICN(JPOS)
+ 130 CONTINUE
+ 140 CONTINUE
+ IF (NZ.GE.NZ1) GO TO 180
+C MOVE UP ENTRIES TO ALLOW FOR DIAGONALS.
+ IPOS = NZ1
+ JPOS = NZ1 - N
+ DO 170 II = 1,N
+ I = N - II + 1
+ J1 = IW2(I)
+ J2 = JPOS
+ IF (J1.GT.JPOS) GO TO 160
+ DO 150 JJ = J1,J2
+ IW(IPOS) = IW(JPOS)
+ A(IPOS) = A(JPOS)
+ IPOS = IPOS - 1
+ JPOS = JPOS - 1
+ 150 CONTINUE
+ 160 IW2(I) = IPOS + 1
+ IPOS = IPOS - 1
+ 170 CONTINUE
+C RUN THROUGH ROWS INSERTING DIAGONAL ENTRIES AND FLAGGING BEGINNING
+C OF EACH ROW BY NEGATING FIRST COLUMN INDEX.
+ 180 DO 190 IOLD = 1,N
+ INEW = PERM(IOLD)
+ JPOS = IW2(INEW) - 1
+ IA = LA - N + IOLD
+ A(JPOS) = A(IA)
+ IW(JPOS) = -IOLD
+ 190 CONTINUE
+C MOVE SORTED MATRIX TO THE END OF THE ARRAYS.
+ IPOS = NZ1
+ IA = LA
+ IIW = LIW
+ DO 200 I = 1,NZ1
+ A(IA) = A(IPOS)
+ IW(IIW) = IW(IPOS)
+ IPOS = IPOS - 1
+ IA = IA - 1
+ IIW = IIW - 1
+ 200 CONTINUE
+ GO TO 230
+C **** ERROR RETURN ****
+ 210 INFO(1) = -3
+ INFO(2) = NZ1
+ GO TO 230
+
+ 220 INFO(1) = -4
+ INFO(2) = NZ1 + N
+C
+ 230 RETURN
+
+ END
+ SUBROUTINE MA27OD(N,NZ,A,LA,IW,LIW,PERM,NSTK,NSTEPS,MAXFRT,NELIM,
+ + IW2,ICNTL,CNTL,INFO)
+C
+C FACTORIZATION SUBROUTINE
+C
+C THIS SUBROUTINE OPERATES ON THE INPUT MATRIX ORDERED BY MA27N/ND AND
+C PRODUCES THE FACTORS OF U AND D ('A'=UTRANSPOSE*D*U) FOR USE IN
+C SUBSEQUENT SOLUTIONS. GAUSSIAN ELIMINATION IS USED WITH PIVOTS
+C CHOSEN FROM THE DIAGONAL. TO ENSURE STABILITY, BLOCK PIVOTS OF
+C ORDER 2 WILL BE USED IF THE DIAGONAL ENTRY IS NOT LARGE ENOUGH.
+C
+C N - MUST BE SET TO THE ORDER OF THE MATRIX. IT IS NOT ALTERED.
+C NZ - MUST BE SET TO THE NUMBER OF NON-ZEROS IN UPPER TRIANGLE OF
+C PERMUTED MATRIX. NOT ALTERED BY MA27O/OD.
+C A - MUST BE SET ON INPUT TO MATRIX HELD BY ROWS REORDERED BY
+C PERMUTATION FROM MA27A/AD IN A(LA-NZ+I),I=1,NZ. ON
+C EXIT FROM MA27O/OD, THE FACTORS OF U AND D ARE HELD IN
+C POSITIONS 1 TO POSFAC-1.
+C LA - LENGTH OF ARRAY A. A VALUE FOR LA
+C SUFFICIENT FOR DEFINITE SYSTEMS
+C WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD.
+C IW - MUST BE SET SO THAT,ON INPUT, IW(LIW-NZ+I),I=1,NZ
+C HOLDS THE COLUMN INDEX OF THE ENTRY IN A(LA-NZ+I). ON EXIT,
+C IW HOLDS INTEGER INDEXING INFORMATION ON THE FACTORS.
+C THE ABSOLUTE VALUE OF THE FIRST ENTRY IN IW WILL BE SET TO
+C THE NUMBER OF BLOCK PIVOTS ACTUALLY USED. THIS MAY BE
+C DIFFERENT FROM NSTEPS SINCE NUMERICAL CONSIDERATIONS
+C MAY PREVENT US CHOOSING A PIVOT AT EACH STAGE. IF THIS ENTRY
+C IN IW IS NEGATIVE, THEN AT LEAST ONE TWO BY TWO
+C PIVOT HAS BEEN USED DURING THE DECOMPOSITION.
+C INTEGER INFORMATION ON EACH BLOCK PIVOT ROW FOLLOWS. FOR
+C EACH BLOCK PIVOT ROW THE COLUMN INDICES ARE PRECEDED BY A
+C COUNT OF THE NUMBER OF ROWS AND COLUMNS IN THE BLOCK PIVOT
+C WHERE, IF ONLY ONE ROW IS PRESENT, ONLY THE NUMBER OF
+C COLUMNS TOGETHER WITH A NEGATIVE FLAG IS HELD. THE FIRST
+C COLUMN INDEX FOR A TWO BY TWO PIVOT IS FLAGGED NEGATIVE.
+C LIW - LENGTH OF ARRAY IW. A VALUE FOR LIW SUFFICIENT FOR
+C DEFINITE SYSTEMS
+C WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD
+C PERM - PERM(I) MUST BE SET TO POSITION OF ROW/COLUMN I IN THE
+C TENTATIVE PIVOT ORDER GENERATED BY MA27A/AD.
+C IT IS NOT ALTERED BY MA27O/OD.
+C NSTK - MUST BE LEFT UNCHANGED SINCE OUTPUT FROM MA27A/AD. NSTK(I)
+C GIVES THE NUMBER OF GENERATED STACK ELEMENTS ASSEMBLED AT
+C STAGE I. IT IS NOT ALTERED BY MA27O/OD.
+C NSTEPS - LENGTH OF ARRAYS NSTK AND NELIM. VALUE IS GIVEN ON OUTPUT
+C FROM MA27A/AD (WILL NEVER EXCEED N). IT IS NOT ALTERED BY
+C MA27O/OD.
+C MAXFRT - NEED NOT BE SET ON INPUT. ON OUTPUT
+C MAXFRT WILL BE SET TO THE MAXIMUM FRONT SIZE ENCOUNTERED
+C DURING THE DECOMPOSITION.
+C NELIM - MUST BE UNCHANGED SINCE OUTPUT FROM MA27A/AD. NELIM(I)
+C GIVES THE NUMBER OF ORIGINAL ROWS ASSEMBLED AT STAGE I.
+C IT IS NOT ALTERED BY MA27O/OD.
+C IW2 - INTEGER ARRAY OF LENGTH N. USED AS WORKSPACE BY MA27O/OD.
+C ALTHOUGH WE COULD HAVE USED A SHORT WORD INTEGER IN THE IBM
+C VERSION, WE HAVE NOT DONE SO BECAUSE THERE IS A SPARE
+C FULL INTEGER ARRAY (USED IN THE SORT BEFORE MA27O/OD)
+C AVAILABLE WHEN MA27O/OD IS CALLED FROM MA27B/BD.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C CNTL is a DOUBLE PRECISION array of length 5, see MA27A/AD.
+C INFO is an INTEGER array of length 20, see MA27A/AD.
+C INFO(1) - INTEGER VARIABLE. DIAGNOSTIC FLAG. A ZERO VALUE ON EXIT
+C INDICATES SUCCESS. POSSIBLE NEGATIVE VALUES ARE ...
+C -3 INSUFFICIENT STORAGE FOR IW.
+C -4 INSUFFICIENT STORAGE FOR A.
+C -5 ZERO PIVOT FOUND IN FACTORIZATION OF DEFINITE MATRIX.
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO,HALF,ONE
+ PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0)
+C ..
+C .. Scalar Arguments ..
+ INTEGER LA,LIW,MAXFRT,N,NSTEPS,NZ
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA),CNTL(5)
+ INTEGER IW(LIW),IW2(N),NELIM(NSTEPS),NSTK(NSTEPS),PERM(N)
+ INTEGER ICNTL(30),INFO(20)
+C ..
+C .. Local Scalars ..
+ DOUBLE PRECISION AMAX,AMULT,AMULT1,AMULT2,DETPIV,RMAX,SWOP,
+ + THRESH,TMAX,UU
+ INTEGER AINPUT,APOS,APOS1,APOS2,APOS3,ASTK,ASTK2,AZERO,I,IASS,
+ + IBEG,IDUMMY,IELL,IEND,IEXCH,IFR,IINPUT,IOLDPS,IORG,IPIV,
+ + IPMNP,IPOS,IROW,ISNPIV,ISTK,ISTK2,ISWOP,IWPOS,IX,IY,J,J1,
+ + J2,JCOL,JDUMMY,JFIRST,JJ,JJ1,JJJ,JLAST,JMAX,JMXMIP,JNEW,
+ + JNEXT,JPIV,JPOS,K,K1,K2,KDUMMY,KK,KMAX,KROW,LAELL,LAPOS2,
+ + LIELL,LNASS,LNPIV,LT,LTOPST,NASS,NBLK,NEWEL,NFRONT,NPIV,
+ + NPIVP1,NTOTPV,NUMASS,NUMORG,NUMSTK,PIVSIZ,POSFAC,POSPV1,
+ + POSPV2
+ INTEGER IDIAG
+C IDIAG IS A TEMPORARY FOR THE DISPLACEMENT FROM THE START OF THE
+C ASSEMBLED MATRIX (OF ORDER IX) OF THE DIAGONAL ENTRY IN ITS ROW IY.
+ INTEGER NTWO,NEIG,NCMPBI,NCMPBR,NRLBDU,NIRBDU
+C ..
+C .. External Subroutines ..
+ EXTERNAL MA27PD
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MAX,MIN
+C ..
+C .. Executable Statements ..
+C INITIALIZATION.
+C NBLK IS THE NUMBER OF BLOCK PIVOTS USED.
+ NBLK = 0
+ NTWO = 0
+ NEIG = 0
+ NCMPBI = 0
+ NCMPBR = 0
+ MAXFRT = 0
+ NRLBDU = 0
+ NIRBDU = 0
+C A PRIVATE VARIABLE UU IS SET TO CNTL(1), SO THAT CNTL(1) WILL REMAIN
+C UNALTERED.
+ UU = MIN(CNTL(1),HALF)
+ UU = MAX(UU,-HALF)
+ DO 10 I = 1,N
+ IW2(I) = 0
+ 10 CONTINUE
+C IWPOS IS POINTER TO FIRST FREE POSITION FOR FACTORS IN IW.
+C POSFAC IS POINTER FOR FACTORS IN A. AT EACH PASS THROUGH THE
+C MAJOR LOOP POSFAC INITIALLY POINTS TO THE FIRST FREE LOCATION
+C IN A AND THEN IS SET TO THE POSITION OF THE CURRENT PIVOT IN A.
+C ISTK IS POINTER TO TOP OF STACK IN IW.
+C ISTK2 IS POINTER TO BOTTOM OF STACK IN IW (NEEDED BY COMPRESS).
+C ASTK IS POINTER TO TOP OF STACK IN A.
+C ASTK2 IS POINTER TO BOTTOM OF STACK IN A (NEEDED BY COMPRESS).
+C IINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN IW.
+C AINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN A.
+C AZERO IS POINTER TO LAST POSITION ZEROED IN A.
+C NTOTPV IS THE TOTAL NUMBER OF PIVOTS SELECTED. THIS IS USED
+C TO DETERMINE WHETHER THE MATRIX IS SINGULAR.
+ IWPOS = 2
+ POSFAC = 1
+ ISTK = LIW - NZ + 1
+ ISTK2 = ISTK - 1
+ ASTK = LA - NZ + 1
+ ASTK2 = ASTK - 1
+ IINPUT = ISTK
+ AINPUT = ASTK
+ AZERO = 0
+ NTOTPV = 0
+C NUMASS IS THE ACCUMULATED NUMBER OF ROWS ASSEMBLED SO FAR.
+ NUMASS = 0
+C
+C EACH PASS THROUGH THIS MAIN LOOP PERFORMS ALL THE OPERATIONS
+C ASSOCIATED WITH ONE SET OF ASSEMBLY/ELIMINATIONS.
+ DO 760 IASS = 1,NSTEPS
+C NASS WILL BE SET TO THE NUMBER OF FULLY ASSEMBLED VARIABLES IN
+C CURRENT NEWLY CREATED ELEMENT.
+ NASS = NELIM(IASS)
+C NEWEL IS A POINTER INTO IW TO CONTROL OUTPUT OF INTEGER INFORMATION
+C FOR NEWLY CREATED ELEMENT.
+ NEWEL = IWPOS + 1
+C SYMBOLICALLY ASSEMBLE INCOMING ROWS AND GENERATED STACK ELEMENTS
+C ORDERING THE RESULTANT ELEMENT ACCORDING TO PERMUTATION PERM. WE
+C ASSEMBLE THE STACK ELEMENTS FIRST BECAUSE THESE WILL ALREADY BE
+C ORDERED.
+C SET HEADER POINTER FOR MERGE OF INDEX LISTS.
+ JFIRST = N + 1
+C INITIALIZE NUMBER OF VARIABLES IN CURRENT FRONT.
+ NFRONT = 0
+ NUMSTK = NSTK(IASS)
+ LTOPST = 1
+ LNASS = 0
+C JUMP IF NO STACK ELEMENTS ARE BEING ASSEMBLED AT THIS STAGE.
+ IF (NUMSTK.EQ.0) GO TO 80
+ J2 = ISTK - 1
+ LNASS = NASS
+ LTOPST = ((IW(ISTK)+1)*IW(ISTK))/2
+ DO 70 IELL = 1,NUMSTK
+C ASSEMBLE ELEMENT IELL PLACING
+C THE INDICES INTO A LINKED LIST IN IW2 ORDERED
+C ACCORDING TO PERM.
+ JNEXT = JFIRST
+ JLAST = N + 1
+ J1 = J2 + 2
+ J2 = J1 - 1 + IW(J1-1)
+C RUN THROUGH INDEX LIST OF STACK ELEMENT IELL.
+ DO 60 JJ = J1,J2
+ J = IW(JJ)
+C JUMP IF ALREADY ASSEMBLED
+ IF (IW2(J).GT.0) GO TO 60
+ JNEW = PERM(J)
+C IF VARIABLE WAS PREVIOUSLY FULLY SUMMED BUT WAS NOT PIVOTED ON
+C EARLIER BECAUSE OF NUMERICAL TEST, INCREMENT NUMBER OF FULLY
+C SUMMED ROWS/COLUMNS IN FRONT.
+ IF (JNEW.LE.NUMASS) NASS = NASS + 1
+C FIND POSITION IN LINKED LIST FOR NEW VARIABLE. NOTE THAT WE START
+C FROM WHERE WE LEFT OFF AFTER ASSEMBLY OF PREVIOUS VARIABLE.
+ DO 20 IDUMMY = 1,N
+ IF (JNEXT.EQ.N+1) GO TO 30
+ IF (PERM(JNEXT).GT.JNEW) GO TO 30
+ JLAST = JNEXT
+ JNEXT = IW2(JLAST)
+ 20 CONTINUE
+ 30 IF (JLAST.NE.N+1) GO TO 40
+ JFIRST = J
+ GO TO 50
+
+ 40 IW2(JLAST) = J
+ 50 IW2(J) = JNEXT
+ JLAST = J
+C INCREMENT NUMBER OF VARIABLES IN THE FRONT.
+ NFRONT = NFRONT + 1
+ 60 CONTINUE
+ 70 CONTINUE
+ LNASS = NASS - LNASS
+C NOW INCORPORATE ORIGINAL ROWS. NOTE THAT THE COLUMNS IN THESE ROWS
+C NEED NOT BE IN ORDER. WE ALSO PERFORM
+C A SWOP SO THAT THE DIAGONAL ENTRY IS THE FIRST IN ITS
+C ROW. THIS ALLOWS US TO AVOID STORING THE INVERSE OF ARRAY PERM.
+ 80 NUMORG = NELIM(IASS)
+ J1 = IINPUT
+ DO 150 IORG = 1,NUMORG
+ J = -IW(J1)
+ DO 140 IDUMMY = 1,LIW
+ JNEW = PERM(J)
+C JUMP IF VARIABLE ALREADY INCLUDED.
+ IF (IW2(J).GT.0) GO TO 130
+C HERE WE MUST ALWAYS START OUR SEARCH AT THE BEGINNING.
+ JLAST = N + 1
+ JNEXT = JFIRST
+ DO 90 JDUMMY = 1,N
+ IF (JNEXT.EQ.N+1) GO TO 100
+ IF (PERM(JNEXT).GT.JNEW) GO TO 100
+ JLAST = JNEXT
+ JNEXT = IW2(JLAST)
+ 90 CONTINUE
+ 100 IF (JLAST.NE.N+1) GO TO 110
+ JFIRST = J
+ GO TO 120
+
+ 110 IW2(JLAST) = J
+ 120 IW2(J) = JNEXT
+C INCREMENT NUMBER OF VARIABLES IN FRONT.
+ NFRONT = NFRONT + 1
+ 130 J1 = J1 + 1
+ IF (J1.GT.LIW) GO TO 150
+ J = IW(J1)
+ IF (J.LT.0) GO TO 150
+ 140 CONTINUE
+ 150 CONTINUE
+C NOW RUN THROUGH LINKED LIST IW2 PUTTING INDICES OF VARIABLES IN NEW
+C ELEMENT INTO IW AND SETTING IW2 ENTRY TO POINT TO THE RELATIVE
+C POSITION OF THE VARIABLE IN THE NEW ELEMENT.
+ IF (NEWEL+NFRONT.LT.ISTK) GO TO 160
+C COMPRESS IW.
+ CALL MA27PD(A,IW,ISTK,ISTK2,IINPUT,2,NCMPBR,NCMPBI)
+ IF (NEWEL+NFRONT.LT.ISTK) GO TO 160
+ INFO(2) = LIW + 1 + NEWEL + NFRONT - ISTK
+ GO TO 770
+
+ 160 J = JFIRST
+ DO 170 IFR = 1,NFRONT
+ NEWEL = NEWEL + 1
+ IW(NEWEL) = J
+ JNEXT = IW2(J)
+ IW2(J) = NEWEL - (IWPOS+1)
+ J = JNEXT
+ 170 CONTINUE
+C
+C ASSEMBLE REALS INTO FRONTAL MATRIX.
+ MAXFRT = MAX(MAXFRT,NFRONT)
+ IW(IWPOS) = NFRONT
+C FIRST ZERO OUT FRONTAL MATRIX AS APPROPRIATE FIRST CHECKING TO SEE
+C IF THERE IS SUFFICIENT SPACE.
+ LAELL = ((NFRONT+1)*NFRONT)/2
+ APOS2 = POSFAC + LAELL - 1
+ IF (NUMSTK.NE.0) LNASS = LNASS* (2*NFRONT-LNASS+1)/2
+ IF (POSFAC+LNASS-1.GE.ASTK) GO TO 180
+ IF (APOS2.LT.ASTK+LTOPST-1) GO TO 190
+C COMPRESS A.
+ 180 CALL MA27PD(A,IW,ASTK,ASTK2,AINPUT,1,NCMPBR,NCMPBI)
+ IF (POSFAC+LNASS-1.GE.ASTK) GO TO 780
+ IF (APOS2.GE.ASTK+LTOPST-1) GO TO 780
+ 190 IF (APOS2.LE.AZERO) GO TO 220
+ APOS = AZERO + 1
+ LAPOS2 = MIN(APOS2,ASTK-1)
+ IF (LAPOS2.LT.APOS) GO TO 210
+ DO 200 K = APOS,LAPOS2
+ A(K) = ZERO
+ 200 CONTINUE
+ 210 AZERO = APOS2
+C JUMP IF THERE ARE NO STACK ELEMENTS TO ASSEMBLE.
+ 220 IF (NUMSTK.EQ.0) GO TO 260
+C PLACE REALS CORRESPONDING TO STACK ELEMENTS IN CORRECT POSITIONS IN A.
+ DO 250 IELL = 1,NUMSTK
+ J1 = ISTK + 1
+ J2 = ISTK + IW(ISTK)
+ DO 240 JJ = J1,J2
+ IROW = IW(JJ)
+ IROW = IW2(IROW)
+ IX = NFRONT
+ IY = IROW
+ IDIAG = ((IY-1)* (2*IX-IY+2))/2
+ APOS = POSFAC + IDIAG
+ DO 230 JJJ = JJ,J2
+ J = IW(JJJ)
+ APOS2 = APOS + IW2(J) - IROW
+ A(APOS2) = A(APOS2) + A(ASTK)
+ A(ASTK) = ZERO
+ ASTK = ASTK + 1
+ 230 CONTINUE
+ 240 CONTINUE
+C INCREMENT STACK POINTER.
+ ISTK = J2 + 1
+ 250 CONTINUE
+C INCORPORATE REALS FROM ORIGINAL ROWS.
+ 260 DO 280 IORG = 1,NUMORG
+ J = -IW(IINPUT)
+C WE CAN DO THIS BECAUSE THE DIAGONAL IS NOW THE FIRST ENTRY.
+ IROW = IW2(J)
+ IX = NFRONT
+ IY = IROW
+ IDIAG = ((IY-1)* (2*IX-IY+2))/2
+ APOS = POSFAC + IDIAG
+C THE FOLLOWING LOOP GOES FROM 1 TO NZ BECAUSE THERE MAY BE DUPLICATES.
+ DO 270 IDUMMY = 1,NZ
+ APOS2 = APOS + IW2(J) - IROW
+ A(APOS2) = A(APOS2) + A(AINPUT)
+ AINPUT = AINPUT + 1
+ IINPUT = IINPUT + 1
+ IF (IINPUT.GT.LIW) GO TO 280
+ J = IW(IINPUT)
+ IF (J.LT.0) GO TO 280
+ 270 CONTINUE
+ 280 CONTINUE
+C RESET IW2 AND NUMASS.
+ NUMASS = NUMASS + NUMORG
+ J1 = IWPOS + 2
+ J2 = IWPOS + NFRONT + 1
+ DO 290 K = J1,J2
+ J = IW(K)
+ IW2(J) = 0
+ 290 CONTINUE
+C PERFORM PIVOTING ON ASSEMBLED ELEMENT.
+C NPIV IS THE NUMBER OF PIVOTS SO FAR SELECTED.
+C LNPIV IS THE NUMBER OF PIVOTS SELECTED AFTER THE LAST PASS THROUGH
+C THE THE FOLLOWING LOOP.
+ LNPIV = -1
+ NPIV = 0
+ DO 650 KDUMMY = 1,NASS
+ IF (NPIV.EQ.NASS) GO TO 660
+ IF (NPIV.EQ.LNPIV) GO TO 660
+ LNPIV = NPIV
+ NPIVP1 = NPIV + 1
+C JPIV IS USED AS A FLAG TO INDICATE WHEN 2 BY 2 PIVOTING HAS OCCURRED
+C SO THAT IPIV IS INCREMENTED CORRECTLY.
+ JPIV = 1
+C NASS IS MAXIMUM POSSIBLE NUMBER OF PIVOTS.
+C WE EITHER TAKE THE DIAGONAL ENTRY OR THE 2 BY 2 PIVOT WITH THE
+C LARGEST OFF-DIAGONAL AT EACH STAGE.
+C EACH PASS THROUGH THIS LOOP TRIES TO CHOOSE ONE PIVOT.
+ DO 640 IPIV = NPIVP1,NASS
+ JPIV = JPIV - 1
+C JUMP IF WE HAVE JUST PROCESSED A 2 BY 2 PIVOT.
+ IF (JPIV.EQ.1) GO TO 640
+ IX = NFRONT-NPIV
+ IY = IPIV-NPIV
+ IDIAG = ((IY-1)* (2*IX-IY+2))/2
+ APOS = POSFAC + IDIAG
+C IF THE USER HAS INDICATED THAT THE MATRIX IS DEFINITE, WE
+C DO NOT NEED TO TEST FOR STABILITY BUT WE DO CHECK TO SEE IF THE
+C PIVOT IS NON-ZERO OR HAS CHANGED SIGN.
+C IF IT IS ZERO, WE EXIT WITH AN ERROR. IF IT HAS CHANGED SIGN
+C AND U WAS SET NEGATIVE, THEN WE AGAIN EXIT IMMEDIATELY. IF THE
+C PIVOT CHANGES SIGN AND U WAS ZERO, WE CONTINUE WITH THE
+C FACTORIZATION BUT PRINT A WARNING MESSAGE ON UNIT ICNTL(2).
+C ISNPIV HOLDS A FLAG FOR THE SIGN OF THE PIVOTS TO DATE SO THAT
+C A SIGN CHANGE WHEN DECOMPOSING AN ALLEGEDLY DEFINITE MATRIX CAN
+C BE DETECTED.
+ IF (UU.GT.ZERO) GO TO 320
+ IF (ABS(A(APOS)).LE.CNTL(3)) GO TO 790
+C JUMP IF THIS IS NOT THE FIRST PIVOT TO BE SELECTED.
+ IF (NTOTPV.GT.0) GO TO 300
+C SET ISNPIV.
+ IF (A(APOS).GT.ZERO) ISNPIV = 1
+ IF (A(APOS).LT.ZERO) ISNPIV = -1
+ 300 IF (A(APOS).GT.ZERO .AND. ISNPIV.EQ.1) GO TO 560
+ IF (A(APOS).LT.ZERO .AND. ISNPIV.EQ.-1) GO TO 560
+ IF (INFO(1).NE.2) INFO(2) = 0
+ INFO(2) = INFO(2) + 1
+ INFO(1) = 2
+ I = NTOTPV + 1
+ IF (ICNTL(2).GT.0 .AND. INFO(2).LE.10) THEN
+ WRITE (ICNTL(2),FMT=310) INFO(1),I
+ END IF
+
+ 310 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
+ + ' *** INFO(1) =',I2,/,' PIVOT',I6,
+ + ' HAS DIFFERENT SIGN FROM THE PREVIOUS ONE')
+
+ ISNPIV = -ISNPIV
+ IF (UU.EQ.ZERO) GO TO 560
+ GO TO 800
+
+ 320 AMAX = ZERO
+ TMAX = AMAX
+C FIND LARGEST ENTRY TO RIGHT OF DIAGONAL IN ROW OF PROSPECTIVE PIVOT
+C IN THE FULLY-SUMMED PART. ALSO RECORD COLUMN OF THIS LARGEST
+C ENTRY.
+ J1 = APOS + 1
+ J2 = APOS + NASS - IPIV
+ IF (J2.LT.J1) GO TO 340
+ DO 330 JJ = J1,J2
+ IF (ABS(A(JJ)).LE.AMAX) GO TO 330
+ JMAX = IPIV + JJ - J1 + 1
+ AMAX = ABS(A(JJ))
+ 330 CONTINUE
+C DO SAME AS ABOVE FOR NON-FULLY-SUMMED PART ONLY HERE WE DO NOT NEED
+C TO RECORD COLUMN SO LOOP IS SIMPLER.
+ 340 J1 = J2 + 1
+ J2 = APOS + NFRONT - IPIV
+ IF (J2.LT.J1) GO TO 360
+ DO 350 JJ = J1,J2
+ TMAX = MAX(ABS(A(JJ)),TMAX)
+ 350 CONTINUE
+C NOW CALCULATE LARGEST ENTRY IN OTHER PART OF ROW.
+ 360 RMAX = MAX(TMAX,AMAX)
+ APOS1 = APOS
+ KK = NFRONT - IPIV
+ LT = IPIV - (NPIV+1)
+ IF (LT.EQ.0) GO TO 380
+ DO 370 K = 1,LT
+ KK = KK + 1
+ APOS1 = APOS1 - KK
+ RMAX = MAX(RMAX,ABS(A(APOS1)))
+ 370 CONTINUE
+C JUMP IF STABILITY TEST SATISFIED.
+ 380 IF (ABS(A(APOS)).GT.MAX(CNTL(3),UU*RMAX)) GO TO 450
+C CHECK BLOCK PIVOT OF ORDER 2 FOR STABILITY.
+ IF (ABS(AMAX).LE.CNTL(3)) GO TO 640
+ IX = NFRONT-NPIV
+ IY = JMAX-NPIV
+ IDIAG = ((IY-1)* (2*IX-IY+2))/2
+ APOS2 = POSFAC + IDIAG
+ DETPIV = A(APOS)*A(APOS2) - AMAX*AMAX
+ THRESH = ABS(DETPIV)
+C SET THRESH TO U TIMES THE RECIPROCAL OF THE MAX-NORM OF THE INVERSE
+C OF THE PROSPECTIVE BLOCK.
+ THRESH = THRESH/ (UU*MAX(ABS(A(APOS))+AMAX,
+ + ABS(A(APOS2))+AMAX))
+C CHECK 2 BY 2 PIVOT FOR STABILITY.
+C FIRST CHECK AGAINST ROW IPIV.
+ IF (THRESH.LE.RMAX) GO TO 640
+C FIND LARGEST ENTRY IN ROW JMAX.
+C FIND MAXIMUM TO THE RIGHT OF THE DIAGONAL.
+ RMAX = ZERO
+ J1 = APOS2 + 1
+ J2 = APOS2 + NFRONT - JMAX
+ IF (J2.LT.J1) GO TO 400
+ DO 390 JJ = J1,J2
+ RMAX = MAX(RMAX,ABS(A(JJ)))
+ 390 CONTINUE
+C NOW CHECK TO THE LEFT OF THE DIAGONAL.
+C WE USE TWO LOOPS TO AVOID TESTING FOR ROW IPIV INSIDE THE LOOP.
+ 400 KK = NFRONT - JMAX + 1
+ APOS3 = APOS2
+ JMXMIP = JMAX - IPIV - 1
+ IF (JMXMIP.EQ.0) GO TO 420
+ DO 410 K = 1,JMXMIP
+ APOS2 = APOS2 - KK
+ KK = KK + 1
+ RMAX = MAX(RMAX,ABS(A(APOS2)))
+ 410 CONTINUE
+ 420 IPMNP = IPIV - NPIV - 1
+ IF (IPMNP.EQ.0) GO TO 440
+ APOS2 = APOS2 - KK
+ KK = KK + 1
+ DO 430 K = 1,IPMNP
+ APOS2 = APOS2 - KK
+ KK = KK + 1
+ RMAX = MAX(RMAX,ABS(A(APOS2)))
+ 430 CONTINUE
+ 440 IF (THRESH.LE.RMAX) GO TO 640
+ PIVSIZ = 2
+ GO TO 460
+
+ 450 PIVSIZ = 1
+ 460 IROW = IPIV - NPIV
+C
+C PIVOT HAS BEEN CHOSEN. IF BLOCK PIVOT OF ORDER 2, PIVSIZ IS EQUAL TO
+C TWO OTHERWISE PIVSIZ EQUALS ONE..
+C THE FOLLOWING LOOP MOVES THE PIVOT BLOCK TO THE TOP LEFT HAND CORNER
+C OF THE FRONTAL MATRIX.
+ DO 550 KROW = 1,PIVSIZ
+C WE JUMP IF SWOP IS NOT NECESSARY.
+ IF (IROW.EQ.1) GO TO 530
+ J1 = POSFAC + IROW
+ J2 = POSFAC + NFRONT - (NPIV+1)
+ IF (J2.LT.J1) GO TO 480
+ APOS2 = APOS + 1
+C SWOP PORTION OF ROWS WHOSE COLUMN INDICES ARE GREATER THAN LATER ROW.
+ DO 470 JJ = J1,J2
+ SWOP = A(APOS2)
+ A(APOS2) = A(JJ)
+ A(JJ) = SWOP
+ APOS2 = APOS2 + 1
+ 470 CONTINUE
+ 480 J1 = POSFAC + 1
+ J2 = POSFAC + IROW - 2
+ APOS2 = APOS
+ KK = NFRONT - (IROW+NPIV)
+ IF (J2.LT.J1) GO TO 500
+C SWOP PORTION OF ROWS/COLUMNS WHOSE INDICES LIE BETWEEN THE TWO ROWS.
+ DO 490 JJJ = J1,J2
+ JJ = J2 - JJJ + J1
+ KK = KK + 1
+ APOS2 = APOS2 - KK
+ SWOP = A(APOS2)
+ A(APOS2) = A(JJ)
+ A(JJ) = SWOP
+ 490 CONTINUE
+ 500 IF (NPIV.EQ.0) GO TO 520
+ APOS1 = POSFAC
+ KK = KK + 1
+ APOS2 = APOS2 - KK
+C SWOP PORTION OF COLUMNS WHOSE INDICES ARE LESS THAN EARLIER ROW.
+ DO 510 JJ = 1,NPIV
+ KK = KK + 1
+ APOS1 = APOS1 - KK
+ APOS2 = APOS2 - KK
+ SWOP = A(APOS2)
+ A(APOS2) = A(APOS1)
+ A(APOS1) = SWOP
+ 510 CONTINUE
+C SWOP DIAGONALS AND INTEGER INDEXING INFORMATION
+ 520 SWOP = A(APOS)
+ A(APOS) = A(POSFAC)
+ A(POSFAC) = SWOP
+ IPOS = IWPOS + NPIV + 2
+ IEXCH = IWPOS + IROW + NPIV + 1
+ ISWOP = IW(IPOS)
+ IW(IPOS) = IW(IEXCH)
+ IW(IEXCH) = ISWOP
+ 530 IF (PIVSIZ.EQ.1) GO TO 550
+C SET VARIABLES FOR THE SWOP OF SECOND ROW OF BLOCK PIVOT.
+ IF (KROW.EQ.2) GO TO 540
+ IROW = JMAX - (NPIV+1)
+ JPOS = POSFAC
+ POSFAC = POSFAC + NFRONT - NPIV
+ NPIV = NPIV + 1
+ APOS = APOS3
+ GO TO 550
+C RESET VARIABLES PREVIOUSLY SET FOR SECOND PASS.
+ 540 NPIV = NPIV - 1
+ POSFAC = JPOS
+ 550 CONTINUE
+C
+ IF (PIVSIZ.EQ.2) GO TO 600
+C PERFORM THE ELIMINATION USING ENTRY (IPIV,IPIV) AS PIVOT.
+C WE STORE U AND DINVERSE.
+ 560 A(POSFAC) = ONE/A(POSFAC)
+ IF (A(POSFAC).LT.ZERO) NEIG = NEIG + 1
+ J1 = POSFAC + 1
+ J2 = POSFAC + NFRONT - (NPIV+1)
+ IF (J2.LT.J1) GO TO 590
+ IBEG = J2 + 1
+ DO 580 JJ = J1,J2
+ AMULT = -A(JJ)*A(POSFAC)
+ IEND = IBEG + NFRONT - (NPIV+JJ-J1+2)
+C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
+CDIR$ IVDEP
+ DO 570 IROW = IBEG,IEND
+ JCOL = JJ + IROW - IBEG
+ A(IROW) = A(IROW) + AMULT*A(JCOL)
+ 570 CONTINUE
+ IBEG = IEND + 1
+ A(JJ) = AMULT
+ 580 CONTINUE
+ 590 NPIV = NPIV + 1
+ NTOTPV = NTOTPV + 1
+ JPIV = 1
+ POSFAC = POSFAC + NFRONT - NPIV + 1
+ GO TO 640
+C PERFORM ELIMINATION USING BLOCK PIVOT OF ORDER TWO.
+C REPLACE BLOCK PIVOT BY ITS INVERSE.
+C SET FLAG TO INDICATE USE OF 2 BY 2 PIVOT IN IW.
+ 600 IPOS = IWPOS + NPIV + 2
+ NTWO = NTWO + 1
+ IW(IPOS) = -IW(IPOS)
+ POSPV1 = POSFAC
+ POSPV2 = POSFAC + NFRONT - NPIV
+ SWOP = A(POSPV2)
+ IF (DETPIV.LT.ZERO) NEIG = NEIG + 1
+ IF (DETPIV.GT.ZERO .AND. SWOP.LT.ZERO) NEIG = NEIG + 2
+ A(POSPV2) = A(POSPV1)/DETPIV
+ A(POSPV1) = SWOP/DETPIV
+ A(POSPV1+1) = -A(POSPV1+1)/DETPIV
+ J1 = POSPV1 + 2
+ J2 = POSPV1 + NFRONT - (NPIV+1)
+ IF (J2.LT.J1) GO TO 630
+ JJ1 = POSPV2
+ IBEG = POSPV2 + NFRONT - (NPIV+1)
+ DO 620 JJ = J1,J2
+ JJ1 = JJ1 + 1
+ AMULT1 = - (A(POSPV1)*A(JJ)+A(POSPV1+1)*A(JJ1))
+ AMULT2 = - (A(POSPV1+1)*A(JJ)+A(POSPV2)*A(JJ1))
+ IEND = IBEG + NFRONT - (NPIV+JJ-J1+3)
+C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
+CDIR$ IVDEP
+ DO 610 IROW = IBEG,IEND
+ K1 = JJ + IROW - IBEG
+ K2 = JJ1 + IROW - IBEG
+ A(IROW) = A(IROW) + AMULT1*A(K1) + AMULT2*A(K2)
+ 610 CONTINUE
+ IBEG = IEND + 1
+ A(JJ) = AMULT1
+ A(JJ1) = AMULT2
+ 620 CONTINUE
+ 630 NPIV = NPIV + 2
+ NTOTPV = NTOTPV + 2
+ JPIV = 2
+ POSFAC = POSPV2 + NFRONT - NPIV + 1
+ 640 CONTINUE
+ 650 CONTINUE
+C END OF MAIN ELIMINATION LOOP.
+C
+ 660 IF (NPIV.NE.0) NBLK = NBLK + 1
+ IOLDPS = IWPOS
+ IWPOS = IWPOS + NFRONT + 2
+ IF (NPIV.EQ.0) GO TO 690
+ IF (NPIV.GT.1) GO TO 680
+ IW(IOLDPS) = -IW(IOLDPS)
+ DO 670 K = 1,NFRONT
+ J1 = IOLDPS + K
+ IW(J1) = IW(J1+1)
+ 670 CONTINUE
+ IWPOS = IWPOS - 1
+ GO TO 690
+
+ 680 IW(IOLDPS+1) = NPIV
+C COPY REMAINDER OF ELEMENT TO TOP OF STACK
+ 690 LIELL = NFRONT - NPIV
+ IF (LIELL.EQ.0 .OR. IASS.EQ.NSTEPS) GO TO 750
+ IF (IWPOS+LIELL.LT.ISTK) GO TO 700
+ CALL MA27PD(A,IW,ISTK,ISTK2,IINPUT,2,NCMPBR,NCMPBI)
+ 700 ISTK = ISTK - LIELL - 1
+ IW(ISTK) = LIELL
+ J1 = ISTK
+ KK = IWPOS - LIELL - 1
+C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
+CDIR$ IVDEP
+ DO 710 K = 1,LIELL
+ J1 = J1 + 1
+ KK = KK + 1
+ IW(J1) = IW(KK)
+ 710 CONTINUE
+C WE COPY IN REVERSE DIRECTION TO AVOID OVERWRITE PROBLEMS.
+ LAELL = ((LIELL+1)*LIELL)/2
+ KK = POSFAC + LAELL
+ IF (KK.NE.ASTK) GO TO 720
+ ASTK = ASTK - LAELL
+ GO TO 740
+C THE MOVE AND ZEROING OF ARRAY A IS PERFORMED WITH TWO LOOPS SO
+C THAT THE CRAY-1 WILL VECTORIZE THEM SAFELY.
+ 720 KMAX = KK - 1
+C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
+CDIR$ IVDEP
+ DO 730 K = 1,LAELL
+ KK = KK - 1
+ ASTK = ASTK - 1
+ A(ASTK) = A(KK)
+ 730 CONTINUE
+ KMAX = MIN(KMAX,ASTK-1)
+ DO 735 K = KK,KMAX
+ A(K) = ZERO
+ 735 CONTINUE
+ 740 AZERO = MIN(AZERO,ASTK-1)
+ 750 IF (NPIV.EQ.0) IWPOS = IOLDPS
+ 760 CONTINUE
+C
+C END OF LOOP ON TREE NODES.
+C
+ IW(1) = NBLK
+ IF (NTWO.GT.0) IW(1) = -NBLK
+ NRLBDU = POSFAC - 1
+ NIRBDU = IWPOS - 1
+ IF (NTOTPV.EQ.N) GO TO 810
+ INFO(1) = 3
+ INFO(2) = NTOTPV
+ GO TO 810
+C **** ERROR RETURNS ****
+ 770 INFO(1) = -3
+ GO TO 810
+
+ 780 INFO(1) = -4
+ INFO(2) = LA + MAX(POSFAC+LNASS,APOS2-LTOPST+2) - ASTK
+ GO TO 810
+
+ 790 INFO(1) = -5
+ INFO(2) = NTOTPV + 1
+ GO TO 810
+
+ 800 INFO(1) = -6
+ INFO(2) = NTOTPV + 1
+ 810 CONTINUE
+ INFO(9) = NRLBDU
+ INFO(10) = NIRBDU
+ INFO(12) = NCMPBR
+ INFO(13) = NCMPBI
+ INFO(14) = NTWO
+ INFO(15) = NEIG
+
+ RETURN
+ END
+ SUBROUTINE MA27PD(A,IW,J1,J2,ITOP,IREAL,NCMPBR,NCMPBI)
+C THIS SUBROUTINE PERFORMS A VERY SIMPLE COMPRESS (BLOCK MOVE).
+C ENTRIES J1 TO J2 (INCL.) IN A OR IW AS APPROPRIATE ARE MOVED TO
+C OCCUPY THE POSITIONS IMMEDIATELY PRIOR TO POSITION ITOP.
+C A/IW HOLD THE ARRAY BEING COMPRESSED.
+C J1/J2 DEFINE THE ENTRIES BEING MOVED.
+C ITOP DEFINES THE POSITION IMMEDIATELY AFTER THE POSITIONS TO WHICH
+C J1 TO J2 ARE MOVED.
+C IREAL MUST BE SET BY THE USER TO 2 IF THE MOVE IS ON ARRAY IW,
+C ANY OTHER VALUE WILL PERFORM THE MOVE ON A.
+C NCMPBR and NCMPBI, see INFO(12) and INFO(13) in MA27A/AD (ACCUMULATE
+C THE NUMBER OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY
+C MA27B/BD.
+C
+C .. Scalar Arguments ..
+ INTEGER IREAL,ITOP,J1,J2,NCMPBR,NCMPBI
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(*)
+ INTEGER IW(*)
+C ..
+C .. Local Scalars ..
+ INTEGER IPOS,JJ,JJJ
+C ..
+C .. Executable Statements ..
+ IPOS = ITOP - 1
+ IF (J2.EQ.IPOS) GO TO 50
+ IF (IREAL.EQ.2) GO TO 20
+ NCMPBR = NCMPBR + 1
+ IF (J1.GT.J2) GO TO 40
+ DO 10 JJJ = J1,J2
+ JJ = J2 - JJJ + J1
+ A(IPOS) = A(JJ)
+ IPOS = IPOS - 1
+ 10 CONTINUE
+ GO TO 40
+
+ 20 NCMPBI = NCMPBI + 1
+ IF (J1.GT.J2) GO TO 40
+ DO 30 JJJ = J1,J2
+ JJ = J2 - JJJ + J1
+ IW(IPOS) = IW(JJ)
+ IPOS = IPOS - 1
+ 30 CONTINUE
+ 40 J2 = ITOP - 1
+ J1 = IPOS + 1
+ 50 RETURN
+
+ END
+ SUBROUTINE MA27QD(N,A,LA,IW,LIW,W,MAXFNT,RHS,IW2,NBLK,LATOP,ICNTL)
+C THIS SUBROUTINE PERFORMS FORWARD ELIMINATION
+C USING THE FACTOR U TRANSPOSE STORED IN A/IA AFTER MA27B/BD.
+C
+C N - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED
+C BY MA27Q/QD.
+C A - MUST BE SET TO HOLD THE REAL VALUES
+C CORRESPONDING TO THE FACTORS OF DINVERSE AND U. THIS MUST BE
+C UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED
+C BY MA27Q/QD.
+C LA - LENGTH OF ARRAY A. NOT ALTERED BY MA27Q/QD.
+C IW - HOLDS THE INTEGER INDEXING
+C INFORMATION FOR THE MATRIX FACTORS IN A. THIS MUST BE
+C UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED
+C BY MA27Q/QD.
+C LIW - LENGTH OF ARRAY IW. NOT ALTERED BY MA27Q/QD.
+C W - USED
+C AS WORKSPACE BY MA27Q/QD TO HOLD THE COMPONENTS OF THE RIGHT
+C HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS.
+C MAXFNT - MUST BE SET TO THE LARGEST NUMBER OF
+C VARIABLES IN ANY BLOCK PIVOT ROW. THIS VALUE WILL HAVE
+C BEEN OUTPUT BY MA27B/BD. NOT ALTERED BY MA27Q/QD.
+C RHS - ON INPUT,
+C MUST BE SET TO HOLD THE RIGHT HAND SIDES FOR THE EQUATIONS
+C WHICH THE USER DESIRES TO SOLVE. ON OUTPUT, RHS WILL HOLD
+C THE MODIFIED VECTORS CORRESPONDING TO PERFORMING FORWARD
+C ELIMINATION ON THE RIGHT HAND SIDES.
+C IW2 - NEED NOT BE SET ON ENTRY. ON EXIT IW2(I) (I = 1,NBLK)
+C WILL HOLD POINTERS TO THE
+C BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW.
+C NBLK - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27Q/QD.
+C LATOP - NEED NOT BE SET ON ENTRY. ON EXIT, IT IS THE POSITION IN
+C A OF THE LAST ENTRY IN THE FACTORS. IT MUST BE PASSED
+C UNCHANGED TO MA27R/RD.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT
+C ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS IS EMPLOYED
+C IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF
+C A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND
+C ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
+C NUMBER OF PIVOTS IN THE BLOCK.
+C
+ INTEGER IFRLVL
+ PARAMETER ( IFRLVL=5 )
+C .. Scalar Arguments ..
+ INTEGER LA,LATOP,LIW,MAXFNT,N,NBLK
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA),RHS(N),W(MAXFNT)
+ INTEGER IW(LIW),IW2(NBLK),ICNTL(30)
+C ..
+C .. Local Scalars ..
+ DOUBLE PRECISION W1,W2
+ INTEGER APOS,IBLK,IFR,ILVL,IPIV,IPOS,IRHS,IROW,IST,J,J1,J2,J3,JJ,
+ + JPIV,K,K1,K2,K3,LIELL,NPIV
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MIN
+C ..
+C .. Executable Statements ..
+C APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A.
+C IPOS. RUNNING POINTER TO BEGINNING OF BLOCK PIVOT ROW IN IW.
+ APOS = 1
+ IPOS = 1
+ J2 = 0
+ IBLK = 0
+ NPIV = 0
+ DO 140 IROW = 1,N
+ IF (NPIV.GT.0) GO TO 90
+ IBLK = IBLK + 1
+ IF (IBLK.GT.NBLK) GO TO 150
+ IPOS = J2 + 1
+C SET UP POINTER IN PREPARATION FOR BACK SUBSTITUTION.
+ IW2(IBLK) = IPOS
+C ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW.
+ LIELL = -IW(IPOS)
+C NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT.
+ NPIV = 1
+ IF (LIELL.GT.0) GO TO 10
+ LIELL = -LIELL
+ IPOS = IPOS + 1
+ NPIV = IW(IPOS)
+ 10 J1 = IPOS + 1
+ J2 = IPOS + LIELL
+ ILVL = MIN(NPIV,10)
+ IF (LIELL.LT.ICNTL(IFRLVL+ILVL)) GO TO 90
+C
+C PERFORM OPERATIONS USING DIRECT ADDRESSING.
+C
+C LOAD APPROPRIATE COMPONENTS OF RIGHT HAND SIDES INTO ARRAY W.
+ IFR = 0
+ DO 20 JJ = J1,J2
+ J = ABS(IW(JJ)+0)
+ IFR = IFR + 1
+ W(IFR) = RHS(J)
+ 20 CONTINUE
+C JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER
+C THE USE OF A 2 BY 2 PIVOT.
+ JPIV = 1
+ J3 = J1
+C PERFORM OPERATIONS.
+ DO 70 IPIV = 1,NPIV
+ JPIV = JPIV - 1
+ IF (JPIV.EQ.1) GO TO 70
+C JUMP IF WE HAVE A 2 BY 2 PIVOT.
+ IF (IW(J3).LT.0) GO TO 40
+C PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT.
+ JPIV = 1
+ J3 = J3 + 1
+ APOS = APOS + 1
+ IST = IPIV + 1
+ IF (LIELL.LT.IST) GO TO 70
+ W1 = W(IPIV)
+ K = APOS
+ DO 30 J = IST,LIELL
+ W(J) = W(J) + A(K)*W1
+ K = K + 1
+ 30 CONTINUE
+ APOS = APOS + LIELL - IST + 1
+ GO TO 70
+C PERFORM OPERATIONS WITH 2 BY 2 PIVOT.
+ 40 JPIV = 2
+ J3 = J3 + 2
+ APOS = APOS + 2
+ IST = IPIV + 2
+ IF (LIELL.LT.IST) GO TO 60
+ W1 = W(IPIV)
+ W2 = W(IPIV+1)
+ K1 = APOS
+ K2 = APOS + LIELL - IPIV
+ DO 50 J = IST,LIELL
+ W(J) = W(J) + W1*A(K1) + W2*A(K2)
+ K1 = K1 + 1
+ K2 = K2 + 1
+ 50 CONTINUE
+ 60 APOS = APOS + 2* (LIELL-IST+1) + 1
+ 70 CONTINUE
+C RELOAD W BACK INTO RHS.
+ IFR = 0
+ DO 80 JJ = J1,J2
+ J = ABS(IW(JJ)+0)
+ IFR = IFR + 1
+ RHS(J) = W(IFR)
+ 80 CONTINUE
+ NPIV = 0
+ GO TO 140
+C
+C PERFORM OPERATIONS USING INDIRECT ADDRESSING.
+C
+C JUMP IF WE HAVE A 2 BY 2 PIVOT.
+ 90 IF (IW(J1).LT.0) GO TO 110
+C PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT.
+ NPIV = NPIV - 1
+ APOS = APOS + 1
+ J1 = J1 + 1
+ IF (J1.GT.J2) GO TO 140
+ IRHS = IW(J1-1)
+ W1 = RHS(IRHS)
+ K = APOS
+ DO 100 J = J1,J2
+ IRHS = ABS(IW(J)+0)
+ RHS(IRHS) = RHS(IRHS) + A(K)*W1
+ K = K + 1
+ 100 CONTINUE
+ APOS = APOS + J2 - J1 + 1
+ GO TO 140
+C PERFORM OPERATIONS WITH 2 BY 2 PIVOT
+ 110 NPIV = NPIV - 2
+ J1 = J1 + 2
+ APOS = APOS + 2
+ IF (J1.GT.J2) GO TO 130
+ IRHS = -IW(J1-2)
+ W1 = RHS(IRHS)
+ IRHS = IW(J1-1)
+ W2 = RHS(IRHS)
+ K1 = APOS
+ K3 = APOS + J2 - J1 + 2
+ DO 120 J = J1,J2
+ IRHS = ABS(IW(J)+0)
+ RHS(IRHS) = RHS(IRHS) + W1*A(K1) + W2*A(K3)
+ K1 = K1 + 1
+ K3 = K3 + 1
+ 120 CONTINUE
+ 130 APOS = APOS + 2* (J2-J1+1) + 1
+ 140 CONTINUE
+ 150 LATOP = APOS - 1
+ RETURN
+
+ END
+ SUBROUTINE MA27RD(N,A,LA,IW,LIW,W,MAXFNT,RHS,IW2,NBLK,LATOP,ICNTL)
+C THIS SUBROUTINE PERFORMS BACKWARD ELIMINATION OPERATIONS
+C USING THE FACTORS DINVERSE AND U
+C STORED IN A/IW AFTER MA27B/BD.
+C
+C N - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED
+C BY MA27R/RD.
+C A - MUST BE SET TO HOLD THE REAL VALUES CORRESPONDING
+C TO THE FACTORS OF DINVERSE AND U. THIS MUST BE
+C UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED
+C BY MA27R/RD.
+C LA - LENGTH OF ARRAY A. NOT ALTERED BY MA27R/RD.
+C IW - HOLDS THE INTEGER INDEXING
+C INFORMATION FOR THE MATRIX FACTORS IN A. THIS MUST BE
+C UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD. NOT ALTERED
+C BY MA27R/RD.
+C LIW - LENGTH OF ARRAY IW. NOT ALTERED BY MA27R/RD.
+C W - USED
+C AS WORKSPACE BY MA27R/RD TO HOLD THE COMPONENTS OF THE RIGHT
+C HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS.
+C MAXFNT - INTEGER VARIABLE. MUST BE SET TO THE LARGEST NUMBER OF
+C VARIABLES IN ANY BLOCK PIVOT ROW. THIS VALUE WAS GIVEN
+C ON OUTPUT FROM MA27B/BD. NOT ALTERED BY MA27R/RD.
+C RHS - ON INPUT,
+C MUST BE SET TO HOLD THE RIGHT HAND SIDE MODIFIED BY THE
+C FORWARD SUBSTITUTION OPERATIONS. ON OUTPUT, RHS WILL HOLD
+C THE SOLUTION VECTOR.
+C IW2 - ON ENTRY IW2(I) (I = 1,NBLK)
+C MUST HOLD POINTERS TO THE
+C BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW, AS SET BY
+C MA27Q/QD.
+C NBLK - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27R/RD.
+C LATOP - IT IS THE POSITION IN
+C A OF THE LAST ENTRY IN THE FACTORS. IT MUST BE UNCHANGED
+C SINCE THE CALL TO MA27Q/QD. IT IS NOT ALTERED BY MA27R/RD.
+C ICNTL is an INTEGER array of length 30, see MA27A/AD.
+C ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT
+C ACCESS IS USED BY MA27C/CD. INDIRECT ACCESS IS EMPLOYED
+C IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF
+C A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND
+C ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
+C NUMBER OF PIVOTS IN THE BLOCK.
+C
+ INTEGER IFRLVL
+ PARAMETER ( IFRLVL=5 )
+C
+C .. Scalar Arguments ..
+ INTEGER LA,LATOP,LIW,MAXFNT,N,NBLK
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LA),RHS(N),W(MAXFNT)
+ INTEGER IW(LIW),IW2(NBLK),ICNTL(30)
+C ..
+C .. Local Scalars ..
+ DOUBLE PRECISION W1,W2
+ INTEGER APOS,APOS2,I1RHS,I2RHS,IBLK,IFR,IIPIV,IIRHS,ILVL,IPIV,
+ + IPOS,IRHS,IST,J,J1,J2,JJ,JJ1,JJ2,JPIV,JPOS,K,LIELL,LOOP,
+ + NPIV
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS,MIN
+C ..
+C .. Executable Statements ..
+C APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A.
+C IPOS. RUNNING POINTER TO BEGINNING OF CURRENT BLOCK PIVOT ROW.
+ APOS = LATOP + 1
+ NPIV = 0
+ IBLK = NBLK + 1
+C RUN THROUGH BLOCK PIVOT ROWS IN THE REVERSE ORDER.
+ DO 180 LOOP = 1,N
+ IF (NPIV.GT.0) GO TO 110
+ IBLK = IBLK - 1
+ IF (IBLK.LT.1) GO TO 190
+ IPOS = IW2(IBLK)
+C ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW.
+ LIELL = -IW(IPOS)
+C NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT.
+ NPIV = 1
+ IF (LIELL.GT.0) GO TO 10
+ LIELL = -LIELL
+ IPOS = IPOS + 1
+ NPIV = IW(IPOS)
+ 10 JPOS = IPOS + NPIV
+ J2 = IPOS + LIELL
+ ILVL = MIN(10,NPIV) + 10
+ IF (LIELL.LT.ICNTL(IFRLVL+ILVL)) GO TO 110
+C
+C PERFORM OPERATIONS USING DIRECT ADDRESSING.
+C
+ J1 = IPOS + 1
+C LOAD APPROPRIATE COMPONENTS OF RHS INTO W.
+ IFR = 0
+ DO 20 JJ = J1,J2
+ J = ABS(IW(JJ)+0)
+ IFR = IFR + 1
+ W(IFR) = RHS(J)
+ 20 CONTINUE
+C JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER
+C THE USE OF A 2 BY 2 PIVOT.
+ JPIV = 1
+C PERFORM ELIMINATIONS.
+ DO 90 IIPIV = 1,NPIV
+ JPIV = JPIV - 1
+ IF (JPIV.EQ.1) GO TO 90
+ IPIV = NPIV - IIPIV + 1
+ IF (IPIV.EQ.1) GO TO 30
+C JUMP IF WE HAVE A 2 BY 2 PIVOT.
+ IF (IW(JPOS-1).LT.0) GO TO 60
+C PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT.
+ 30 JPIV = 1
+ APOS = APOS - (LIELL+1-IPIV)
+ IST = IPIV + 1
+ W1 = W(IPIV)*A(APOS)
+ IF (LIELL.LT.IST) GO TO 50
+ JJ1 = APOS + 1
+ DO 40 J = IST,LIELL
+ W1 = W1 + A(JJ1)*W(J)
+ JJ1 = JJ1 + 1
+ 40 CONTINUE
+ 50 W(IPIV) = W1
+ JPOS = JPOS - 1
+ GO TO 90
+C PERFORM BACK-SUBSTITUTION OPERATIONS WITH 2 BY 2 PIVOT
+ 60 JPIV = 2
+ APOS2 = APOS - (LIELL+1-IPIV)
+ APOS = APOS2 - (LIELL+2-IPIV)
+ IST = IPIV + 1
+ W1 = W(IPIV-1)*A(APOS) + W(IPIV)*A(APOS+1)
+ W2 = W(IPIV-1)*A(APOS+1) + W(IPIV)*A(APOS2)
+ IF (LIELL.LT.IST) GO TO 80
+ JJ1 = APOS + 2
+ JJ2 = APOS2 + 1
+ DO 70 J = IST,LIELL
+ W1 = W1 + W(J)*A(JJ1)
+ W2 = W2 + W(J)*A(JJ2)
+ JJ1 = JJ1 + 1
+ JJ2 = JJ2 + 1
+ 70 CONTINUE
+ 80 W(IPIV-1) = W1
+ W(IPIV) = W2
+ JPOS = JPOS - 2
+ 90 CONTINUE
+C RELOAD WORKING VECTOR INTO SOLUTION VECTOR.
+ IFR = 0
+ DO 100 JJ = J1,J2
+ J = ABS(IW(JJ)+0)
+ IFR = IFR + 1
+ RHS(J) = W(IFR)
+ 100 CONTINUE
+ NPIV = 0
+ GO TO 180
+C
+C PERFORM OPERATIONS USING INDIRECT ADDRESSING.
+C
+ 110 IF (NPIV.EQ.1) GO TO 120
+C JUMP IF WE HAVE A 2 BY 2 PIVOT.
+ IF (IW(JPOS-1).LT.0) GO TO 150
+C PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT.
+ 120 NPIV = NPIV - 1
+ APOS = APOS - (J2-JPOS+1)
+ IIRHS = IW(JPOS)
+ W1 = RHS(IIRHS)*A(APOS)
+ J1 = JPOS + 1
+ IF (J1.GT.J2) GO TO 140
+ K = APOS + 1
+ DO 130 J = J1,J2
+ IRHS = ABS(IW(J)+0)
+ W1 = W1 + A(K)*RHS(IRHS)
+ K = K + 1
+ 130 CONTINUE
+ 140 RHS(IIRHS) = W1
+ JPOS = JPOS - 1
+ GO TO 180
+C PERFORM OPERATIONS WITH 2 BY 2 PIVOT
+ 150 NPIV = NPIV - 2
+ APOS2 = APOS - (J2-JPOS+1)
+ APOS = APOS2 - (J2-JPOS+2)
+ I1RHS = -IW(JPOS-1)
+ I2RHS = IW(JPOS)
+ W1 = RHS(I1RHS)*A(APOS) + RHS(I2RHS)*A(APOS+1)
+ W2 = RHS(I1RHS)*A(APOS+1) + RHS(I2RHS)*A(APOS2)
+ J1 = JPOS + 1
+ IF (J1.GT.J2) GO TO 170
+ JJ1 = APOS + 2
+ JJ2 = APOS2 + 1
+ DO 160 J = J1,J2
+ IRHS = ABS(IW(J)+0)
+ W1 = W1 + RHS(IRHS)*A(JJ1)
+ W2 = W2 + RHS(IRHS)*A(JJ2)
+ JJ1 = JJ1 + 1
+ JJ2 = JJ2 + 1
+ 160 CONTINUE
+ 170 RHS(I1RHS) = W1
+ RHS(I2RHS) = W2
+ JPOS = JPOS - 2
+ 180 CONTINUE
+ 190 RETURN
+
+ END
diff --git a/docker/requirements2.txt b/docker/requirements2.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2aec83576f2d0d46825646b9a59926e3c874ff6a
--- /dev/null
+++ b/docker/requirements2.txt
@@ -0,0 +1,4 @@
+-e git+https://github.com/FUSED-Wind/windIO.git#egg=windio
+-e git+https://github.com/OpenMDAO/OpenMDAO.git#egg=openmdao
+bokeh==0.12.14
+ipywidgets