MINICERN/packlib/kernlib/kerngen/tcgen/tlerr.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/02/15 17:49:53  mclareni
   6 * Kernlib
   7 *
   8 *
   9 #include "kerngen/pilot.h"
  10       SUBROUTINE TLERR (A,X,AUX,IPIV)
  11 C
  12 C CERN PROGLIB# E230    TLERR           .VERSION KERNFOR  2.06  740511
  13 C ORIG. 11/05/74 WH+WM
  14 C
  15 C.  SUBROUTINE TLERR          L.S. ERROR MATRIX             HART/MATT
  16 C.
  17 C.        CALCULATES INVERSE OF (A-TRANSPOSED*A) DIRECTLY FROM THE
  18 C.        TRIANGULARISED TRANSFORM OF A.
  19 C.  ARGUMENTS
  20 C.        A,X,AUX,IPIV,(M1,M,N,L,IER) DEFINED AS FOR TLSC WITH X THE
  21 C.        COVARIANCE MATRIX. MATRIX X MAY OVERWRITE MATRIX A.
  22 C.  REMARKS
  23 C.        CONSTRAINED AND UNCONSTRAINED VERSION COMBINED
  24 C.
  25 C.-------------------------------------------------------------------
  26 C
  27       COMMON /TLSDIM/ M1,M,N,L,IER
  28       COMMON /SLATE/ BETA,H,I,IA,ID,IEND,II,IK,IL,IST,J,JA,JD,JK,JL
  29      +              ,K,KN,KS,K1,LV,NK,PIV,SIG,DUM(17)
  30       DIMENSION      A(*), AUX(*), IPIV(*), X(*)
  31 C
  32       K1 = MAX (N,L)
  33       IF     (IABS (IER).EQ.N)         GO TO     5
  34 C--
  35 C--      COMPLETE HOUSEHOLDER TRANSFORMATION IF IER LESS THAN N.
  36 C
  37       IST = IER * (N+1) + 1
  38       KS  = IER + 1
  39 C
  40       DO           4         K=KS,N
  41       LV = M - K + 1
  42 C
  43 C--      GENERATE VECTOR UK AND TRANSFORMATION PARAMETER BETA.
  44 C
  45       CALL TLUK (A(IST),N,LV,SIG,BETA)
  46       J = K1 + K
  47       AUX(J) = -SIG
  48       IF     (K.EQ.N)        GO TO     4
  49 C
  50 C--      TRANSFORMATION OF MATRIX A.
  51 C--
  52       NK = N - K
  53       IF (LV.EQ.1)                     GO TO     2
  54       CALL TLSTEP (A(IST),A(IST+1),N,N,LV,NK,BETA)
  55       GO TO        4
  56     2 DO           3         J=1,NK
  57       JST = IST + J
  58     3 A(JST) = A(JST)*(1.-BETA*A(IST)**2)
  59       IST = IST + N + 1
  60     4 IPIV(K) = K
  61 C
  62 C
  63 C--      COMPUTE X FROM A AND DIAGONAL ELEMENTS OF A-TRANSPOSED.
  64 C
  65     5 DO           40        K=1,N
  66       KN = N-K+1
  67       IA = (KN-1)*( N+1)+1
  68       IK = KN *N
  69       IL = N*N - K + 1
  70       II = KN + K1
  71       PIV=1./AUX(II)
  72       ID = IPIV(KN)-KN
  73       JA = IA+1
  74       JK = IK
  75       JL = IL
  76 C
  77       DO           20        J=1,K
  78       H=0.
  79       IF     (J.EQ.K .AND. J.LE.N-M1)            H = PIV
  80       IF     (K.EQ.1)                            GO TO     15
  81       II = JK
  82 C
  83       DO           10        I=JA,IK
  84       II=II+N
  85    10 H = H-A(I)*X(II)
  86 C
  87    15 H = H*PIV
  88       X(JL) = H
  89       JK = JK - 1
  90    20 JL = JL - N
  91 C
  92 C--      COMPLETE SYMMETRIC PART.
  93 C
  94       IF     (K.EQ.1)                  GO TO     40
  95       JL = IA
  96       DO           25        J=JA,IK
  97       JL = JL + N
  98    25 X(J) = X(JL)
  99 C
 100 C--      INTERCHANGE OF ROWS AND COLUMNS ALREADY FINISHED.
 101 C
 102       IF     (ID.EQ.0)                 GO TO     40
 103       DO           30        J=IA,IL,N
 104       II = J + ID
 105       H  = X(II)
 106       X(II) = X(J)
 107    30 X(J)  = H
 108 C
 109       ID = ID*N
 110       DO           35        J=IA,IK
 111       II = J + ID
 112       H  = X(II)
 113       X(II) = X(J)
 114    35 X(J)  = H
 115 C
 116    40 CONTINUE
 117 C
 118       RETURN
 119       END