C @(#)cracow.for	17.1.1.1 (ES0-DMD) 01/25/02 17:20:20
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  @(#)cracow.for	5.1 (ESO-IPG) 4/5/93 15:58:27
C
C
C
C   C R A C O W
C
C   Purpose:
C   To find least squares solution of normal equations.
C
C   Input:
C   A(n1,n1) - a symmetric matrix. If matrix A represents normal
C     equations, put the sum of squares of RHS into element A(n1,n1)
C     for calculation of the chi2 and covariance matrix. Only A(i,j)
C     for i>=j are requires for input.
C   N1       - size of matrix A (number of unknowns) +1
C   M        - M>=N1 - number of condition equations
C            - M< N1 - only triangle root and its triangle inverse are
C              returned
C
C   Output:
C   M>=N1:
C     A(i,n1) - unknowns (i=1,...,n)
C     A(i,i)  - their variances
C     A(n1,n1)- variance of an observation of unit weight (Chi2/D.F.)
C     A(i,j)  - for j>=i half of the covariance matrix of unknowns
C
C   M<N1:
C     A(i,j)  - for j>=i inverse triangle root matrix (its square is A**(-1))
C     A(i,j)  - for j<i contains square root matrix except of diagonal.
C               The diagonal contains arithmetic inverses of its elements.
C
C   Method:
C   Banachiewicz-Cholesky decomposition of the normal equations.
C   This subroutine is particularly suitable for well posed problems
C   with near-orthogonal base functions.
C     This method should not be used for near singular matrices,
C   since it sufferes from large numerical errors then. Use the singular
C   value decomposition method in such cases.
C
C   Alex Schwarzenberg-Czerny                         Warsaw, October 1987
C
      DOUBLE PRECISION FUNCTION CRACOW(A,N1,MAXN1,M)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
C
      INTEGER N1,MAXN1,M
      DOUBLE PRECISION A(MAXN1,MAXN1)
C
      INTEGER I,J,K
      DOUBLE PRECISION AII,S
C
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
C
C   Find Banachiewicz-Cholesky triangle root and its inverse Q, so that
C   Q*Q=A**(-1). The unknowns -a(i,n1) are a byproduct. So...
C
      CRACOW=ZERO
      IF (N1.GT.MAXN1) GOTO 99
      DO 1 I=1,N1-1
        IF (A(I,I).LE.ZERO) GOTO 99
        AII=SQRT(A(I,I))
C
C   ... complete previous i step, ...
C
        DO 5 J=1,N1
          A(J,I)=A(J,I)/AII
 5      CONTINUE
C
C   ... complete next row of the triangle matrix, ...
C
        DO 4 J=I+1,N1
          DO 41 K=1,I
            A(J,I+1)=A(J,I+1)-A(J,K)*A(I+1,K)
 41       CONTINUE
 4      CONTINUE
C
C   ... compute next column of the inverse matrix.
C
        A(I,I)=ONE/AII
        DO 6 J=1,I
          A(J,I+1)=ZERO
          DO 7 K=J,I
            A(J,I+1)=A(J,I+1)-A(J,K)*A(I+1,K)
 7        CONTINUE
 6      CONTINUE
 1    CONTINUE
C
C   Compute covariance matrix and reverse sign of unknowns
C
      IF (M.GE.N1) THEN
        A(N1,N1)=A(N1,N1)/(M-N1+1)
        DO 10 I=1,N1-1
          A(I,N1)=-A(I,N1)
          DO 12 J=1,I
            S=ZERO
            DO 11 K=I,N1-1
              S=S+A(I,K)*A(J,K)
 11         CONTINUE
            A(J,I)=S*A(N1,N1)
 12       CONTINUE
 10     CONTINUE
      ENDIF
      CRACOW=ONE
 99   END
C
C
C