MINICERN/mathlib/gen/d/mconv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:21  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE MCONV (N)
  11 C-----MCONV INVERTS THE POSITIVE DEFINITE PACKED SYMMETRIC MATRIX Z
  12 C-----BY THE SQUARE-ROOT METHOD
  13 #include "d510pl.inc"
  14 #include "d510si.inc"
  15 #include "d510uo.inc"
  16 C-----MAXIMUM REAL NUMBER AND 10.*MAXIMUM RELATIVE PRECISION ON CDC6000
  17 #if defined(CERNLIB_CDC)
  18       DATA AM,RP/1.E300,1.E-14/
  19 #endif
  20 #if defined(CERNLIB_IBM)
  21       DATA AM,RP/1.E75,1.E-14/
  22 #endif
  23 #if (!defined(CERNLIB_CDC))&&(!defined(CERNLIB_IBM))
  24       DATA AM,RP / 1.0E37, 1.0E-14/
  25 #endif
  26       IF (N.LT.1) RETURN
  27       APS=SQRT(AM/N)
  28       AP=1./(APS*APS)
  29       IR=0
  30       DO 11 I=1,N
  31  1    IR=IR+1
  32       IF (PL(IR)) 1,1,2
  33  2    NI=I*(I-1)/2
  34       II=NI+I
  35       K=N+1
  36       IF (Z(II).LE.RP*ABS(R(IR)).OR.Z(II).LE.AP) GO TO 19
  37       Z(II)=1./SQRT(Z(II))
  38       NL=II-1
  39  3    IF (NL-NI) 5,5,4
  40  4    Z(NL)=Z(NL)*Z(II)
  41       IF (ABS(Z(NL)).GE.APS) GO TO 16
  42       NL=NL-1
  43       GO TO 3
  44  5    IF (I-N) 6,12,12
  45  6    K=K-1
  46       NK=K*(K-1)/2
  47       NL=NK
  48       KK=NK+I
  49       D=Z(KK)*Z(II)
  50       C=D*Z(II)
  51       L=K
  52  7    LL=NK+L
  53       LI=NL+I
  54       Z(LL)=Z(LL)-Z(LI)*C
  55       L=L-1
  56       NL=NL-L
  57       IF (L-I) 9,9,7
  58  8    LL=NK+L
  59       LI=NI+L
  60       Z(LL)=Z(LL)-Z(LI)*D
  61  9    L=L-1
  62       IF (L) 10,10,8
  63  10   Z(KK)=-C
  64       IF (K-I-1) 11,11,6
  65  11   CONTINUE
  66  12   DO 14 I=1,N
  67       DO 14 K=I,N
  68       NL=K*(K-1)/2
  69       KI=NL+I
  70       D=0.
  71       DO 13 L=K,N
  72       LI=NL+I
  73       LK=NL+K
  74       D=D+Z(LI)*Z(LK)
  75       NL=NL+L
  76  13   CONTINUE
  77       KI=K*(K-1)/2+I
  78       Z(KI)=D
  79  14   CONTINUE
  80  15   RETURN
  81  16   K=I+NL-II
  82       IR=0
  83       DO 18 I=1,K
  84  17   IR=IR+1
  85       IF (PL(IR)) 17,17,18
  86  18   CONTINUE
  87  19   PL(IR)=-2.
  88       R(IR)=0.
  89       INDFLG(1)=IR
  90       GO TO 15
  91       END