C @(#)partit.for	17.1.1.1 (ES0-DMD) 01/25/02 17:17:17
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C  @(#)partit.for	17.1.1.1  (ESO-IPG)  01/25/02  17:17:17
      SUBROUTINE PARTIT(A,N,X,NBLOKS,MBLOK,*)
C
C       Copyright (C) Andrew T. Young, 1990, 1991
C
C
C  Inverts compact sym.pos.def.matrix  A (N X N) by partitioning.
C  See Frazer, Duncan & Collar "Elementary Matrices" p.113
C
C  Assume upper left  NP = (NBLOKS * MBLOK) is block-diagonal,
C  with  (NBLOKS)  blocks, each  (MBLOK X MBLOK).
C
C  :     U        V     : -1      :  -1       -1            -1 :
C  :  (P X P)  (P X Q)  :         : U   + X TH  Y     - X TH   :
C  :                    :    =    :                            :
C  :      T             :         :       -1             -1    :
C  :     V        W     :         :   - TH  Y          TH      :
C  :  (Q X P)  (Q X Q)  :         :                            :
C
C  WHERE      -1            T -1       T
C        X = U  V ,    Y = V U   ,   (X  = Y)
C
C                                    T T
C       TH = W - YV = W - S,    S = V Y
C
C                -1      T   -1
C        P = X(TH  Y) = Y (TH  Y)
C
C
c     IMPLICIT NONE
      IMPLICIT INTEGER(I-N)
C
      DOUBLE PRECISION A,X
      INTEGER N,NBLOKS,MBLOK
C
      DIMENSION A(N*(N+1)/2), X(N)
      CHARACTER*60 CARD
C
C   X IS SCRATCH VECTOR.
C
      INTEGER M1
      COMMON/ROWS/ M1
C
      INTEGER NP,NPP1,LVTZ,NROW, J, KR,
     1    NBLOK, MAXK, LRZ, LWR, KROW, K, NRUZ, KOL, LY, NROWU, LYZ, LW 
C
C
      NP=NBLOKS*MBLOK
      NPP1=NP+1
      LVTZ=NP*(NPP1)/2
C
C  COMPUTE U-1, Y, AND THETA.
C
      NROW=1
      DO 20 NBLOK=1,NBLOKS
C   INVERT BLOCK.
      CALL BLKINV(A,MBLOK,NROW,X,*999)
      IF(NP.EQ.N) GO TO 20
C   M1, NROWZ NOW SET FOR THIS BLOCK.
C
C  COMPUTE VT(U-1).
C
C  METHOD:
C   SAVE COLS.(NROW-MAXK) OF VT IN X.  CLEAR SPACE IN VT.
C   COL. KOL OF THIS ROW OF VT(U-1) IS DOT PRODUCT
C   OF THIS ROW (SAVED IN X) WITH COL. KOL OF (U-1),
C   WHICH IS BENT BY TRIANGULAR STORAGE.  1ST PART
C   OF THIS COL. IS IN ROW (KOL-NROW) OF A.
C
      MAXK=NROW+MBLOK-1
      LRZ=LVTZ
      LWR=LVTZ+NPP1
      DO 10 KROW=NPP1,N
C
C   SAVE PART OF NEXT ROW OF VT (KROW OF A) IN X.
      DO 2 K=NROW,MAXK
      X(K)=A(LRZ+K)
    2 A(LRZ+K)=0.
C
      NRUZ=M1-NROW
      DO 7 KOL=NROW,MAXK
      LY=LRZ+KOL
C   COMPUTE 1ST PART OF (KROW,KOL) ELEMENT OF (VT.U-1) VIA ROW OF (U-1).
      DO 4 K=NROW,KOL
    4 A(LY)=A(LY) +  X(K) * A(NRUZ+K)
C   2ND PART OF ROW.  USE COL.OF (U-1).
      NROWU=NRUZ+KOL
      DO 5 K=KOL+1,MAXK
      A(LY)=A(LY) +  X(K) * A(NROWU+KOL)
C     Y(KROW,KOL)   VT(KROW,K)*(U-1)(K,KOL)
    5 NROWU=NROWU+K
    7 NRUZ=NRUZ+KOL
C
C  NOW GET S = (VT)(YT).
C
C   WE STILL HAVE COLS.NCOL TO MAXK OF VT IN X.
C   DOT IT INTO ROW J OF Y TO GET SOME OF S(KROW,J).
C   N.B.:  J .LE. KROW.    EACH BLOCK ADDS TO ALL OF W.
      LYZ=LVTZ
      LW=LWR
      DO 9 J=NPP1,KROW
C   LOOP OVER J = ROWS OF Y = COLS.OF S.
      DO 8 KOL=NROW,MAXK
    8 A(LW)=A(LW)-X(KOL)*A(LYZ+KOL)
C     THETA = W  - (VT) * (YT)
      LYZ=LYZ+J
    9 LW=LW+1
      LWR=LWR+KROW
   10 LRZ=LRZ+KROW
C    THIS ENDS BLOCK NBLOK.
   20 NROW=NROW+MBLOK
      IF(NP.EQ.N) RETURN
C
C   WE NOW HAVE Y AND THETA.  INVERT THETA.
C
      CALL BLKINV(A,N-NP,NPP1,X,*998)
      IF(NP.EQ.0) RETURN
C
C  COMPUTE -(THETA-1)*Y AND P.
C
C   DO COLUMNS.
      DO 30 KOL=1,NP
C   SAVE COL. KOL OF Y IN X.
      LRZ=LVTZ
      DO 21 KR=NPP1,N
      X(KR)=A(LRZ+KOL)
      A(LRZ+KOL)=0.
   21 LRZ=LRZ+KR
C
C   COMPUTE COL.OF (THETA-1)*Y.
      LRZ=LVTZ
      DO 25 KROW=NPP1,N
C   (KROW,KOL) ELEMENT.
      LY=LRZ+KOL
C   1ST PART.
      DO 22 KR=NPP1,KROW
   22 A(LY)=A(LY) - X(KR)*A(LRZ+KR)
C      P(KR,KOL)      Y * (THETA-1)
      LWR=LRZ+KROW+KROW
C   2ND PART.
      DO 23 KR=KROW+1,N
      A(LY)=A(LY) - X(KR)*A(LWR)
C      P(KR,KOL)      Y * (THETA-1)
   23 LWR=LWR+KR
   25 LRZ=LRZ+KROW
C  - (THETA-1)*Y NOW DONE.
C
C  COMPUTE -P = (YT)*((THETA-1)*Y) AND SUBTRACT FROM (U-1).
C
C   X-ARRAY STILL HAS KOL OF Y.  DOT INTO COL.KOL OF -(THETA-1)*Y
C   TO GET -P(KOL,KOL).  USE COL.K (.GT.KOL) OF Y TO GET ROW K OF -P.
C
      LY=KOL*(KOL+1)/2
      LRZ=LVTZ
      DO 27 J=NPP1,N
      A(LY)=A(LY) - X(J)*A(LRZ+KOL)
   27 LRZ=LRZ+J
      LY=LY+KOL
      DO 29 K=KOL+1,NP
      LRZ=LVTZ
      DO 28 J=NPP1,N
      A(LY)=A(LY) - A(LRZ+K)*A(LRZ+KOL)
C    (U-1)(K,KOL)    Y(J,K) * -(THETA-1)*Y (J,KOL)
   28 LRZ=LRZ+J
   29 LY=LY+K
   30 CONTINUE
      RETURN
C
  998 CALL TV('Final block singular:')
      NBLOK=NBLOKS+1
  999 WRITE(CARD,'(A,I3,A,I3,A)')
     1   'PARTIT failed for nblok = ',NBLOK,' of ',NP,' parameters'
      CALL TV(CARD)
      RETURN 1
      END
      SUBROUTINE BLKINV(A,N,NROW1,X,*)
C
C  Alg.9 of Nash, to invert a triang.pos.def.sym.matrix in situ, p.84.
C  modified to invert submatrix of size N starting in row  NROW1.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A,X
      INTEGER N, NROW1
C
      DIMENSION A((N+NROW1)*(N+NROW1-1)/2), X(N)
C
C   X IS SCRATCH SPACE.
C
      INTEGER M1
      COMMON/ROWS/ M1
C
      DOUBLE PRECISION S,T
      INTEGER NROWZ, K, M, I, LQ, LP, J
C
C
      M1=NROW1*(NROW1+1)/2
      NROWZ=NROW1-1
      DO 12 K=N,1,-1
      S=A(M1)
      IF(S.LE.0.) RETURN 1
      M=M1
      DO 9 I=2,N
      LQ=M+NROWZ
      M=LQ+I
C   M IS END OF ROW I; LQ IS JUST BEFORE ROW I.
      T=A(LQ+1)
      X(I)=-T/S
      IF(I.GT.K) X(I)=-X(I)
      LP=NROWZ+I
      DO 8 J=LQ+2,M
    8 A(J-LP)=A(J)+T*X(J-LQ)
    9 CONTINUE
      LQ=LQ-1
      A(M)=1./S
      DO 11 I=2,N
   11 A(LQ+I)=X(I)
   12 CONTINUE
      RETURN
      END