ma27ad.f

* 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)