MINICERN/mathlib/gen/d/linsq.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 LINSQ(K,N,M,F,X,Y,W,DA,H,COV,QPRINT,QMIN)
  11       DIMENSION X(M,N),F(K),Y(N),W(N),H(K,K),DA(K)
  12       DIMENSION INDEX(100)
  13       DOUBLE PRECISION H,DETERM
  14       DOUBLE PRECISION DWI,DYI,DFK1
  15       KPLUS=K+1
  16       APRINT=QPRINT
  17       DO 11 I=1,K
  18       DO 11 J=1,KPLUS
  19    11 H(I,J) = 0.D0
  20       IFLAG=1
  21       DO 2 I=1,N
  22       CALL FCN(K,M,F,X(1,I),IFLAG)
  23       IFLAG=4
  24       DWI = DBLE(W(I))
  25       DYI = DBLE(Y(I))
  26       DO 2 K1=1,K
  27       DFK1 = DBLE(F(K1))
  28       H(K1,K+1) = H(K1,K+1) + DWI*DFK1*DYI
  29       DO 2 L=1,K1
  30     2 H(K1,L) = H(K1,L) + DWI*DFK1*DBLE(F(L))
  31       DO 3 J=1,K
  32       DO 3 I=1,J
  33     3 H(I,J)=H(J,I)
  34       IF(COV)5,4,5
  35     4 CALL D508R2(H,K,K,K+1,1,INDEX,NERROR,DETERM)
  36       MM=1
  37       QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
  38       DO 7 I=1,K
  39     7 F(I)=H(I,1)
  40       IF(APRINT.NE.0.) GO TO 10
  41       RETURN
  42     5 CALL D508R1(H,K,K,K+1,1,INDEX,NERROR,DETERM)
  43       MM=K+1
  44       QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
  45       DO 6 I=1,K
  46     6 F(I)=H(I,K+1)
  47       DO 8 I=1,K
  48       HDA=H(I,I)
  49     8 DA(I)=SQRT(HDA)
  50       IF(APRINT.EQ.0.) RETURN
  51 C
  52 C     PRINTS ERRORS IN THE COEFFICIENTS
  53 C
  54       WRITE(6,100)
  55   100 FORMAT(45X,'ERRORS IN THE COEFFICIENTS')
  56       WRITE(6,101)(DA(I),I=1,K)
  57   101 FORMAT(5X,5G16.6)
  58 C
  59 C     PRINTS LOWER DIAGONAL COVARIANCE MATRIX
  60 C
  61       WRITE(6,102)
  62   102 FORMAT(///40X,'VARIANCE-COVARIANCE MATRIX'///)
  63       DO 9 I=1,K
  64       DO 12 J=1,K
  65    12 F(J)=H(I,J)
  66     9 WRITE(6,103)(F(JJ),JJ=1,I)
  67   103 FORMAT(10X,5G20.6//)
  68       DO 13 I=1,K
  69    13 F(I)=H(I,K+1)
  70    10 CONTINUE
  71 C
  72 C     PRINTS COEFFICIENTS AND SUM OF SQUARES
  73 C
  74       WRITE(6,104)QMIN
  75   104 FORMAT(///30X,'SUM OF SQUARES=',G20.6///)
  76       WRITE(6,105)
  77   105 FORMAT(45X,'COEFFICIENTS'///)
  78       WRITE(6,103)(F(I),I=1,K)
  79       RETURN
  80       END
  81       FUNCTION QLINSQ(K,N,M,F,X,Y,W,DA,H,MM)
  82 C     COMPUTES THE SUM OF SQUARES
  83       DIMENSION F(K),X(M,N),DA(K),H(K,K),Y(N),W(N)
  84       DOUBLE PRECISION H
  85       IFLAG = 4
  86       Q=0.
  87       DO 3 I=1,N
  88       CALL FCN(K,M,F,X(1,I),IFLAG)
  89       HSUM=0.
  90       DO 2 I1=1,K
  91     2 HSUM=H(I1,MM)*F(I1)+HSUM
  92     3 Q=(HSUM-Y(I))*(HSUM-Y(I))*W(I)+Q
  93       QLINSQ = Q
  94       RETURN
  95       END
  96       SUBROUTINE  D508R1 (A,IDIM1,N1,IDIM2,N2,INDEX,NERROR,DETERM)
  97 C
  98 C        MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS
  99       DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
 100       DIMENSION A(IDIM1),INDEX(IDIM1)
 101       DETER=1.0D0
 102       N=N1
 103       IEMAT=N+N2
 104       IDIM=IDIM1
 105       NMIN1=N-1
 106 C        THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF
 107 C        ARRAY A.
 108       IPIVC=1-IDIM
 109 C        MAIN LOOP TO INVERT THE MATRIX
 110       DO 11 MAIN=1,N
 111       PIVOT=0.0D0
 112       IPIVC=IPIVC+IDIM
 113 C        SEARCH FOR NEXT PIVOT IN COLUMN MAIN.
 114       IPIVC1=IPIVC+MAIN-1
 115       IPIVC2=IPIVC +NMIN1
 116       DO 2 I1=IPIVC1,IPIVC2
 117       IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
 118     1 PIVOT=A(I1)
 119       LPIV=I1
 120     2 CONTINUE
 121 C        IS PIVOT DIFFERENT FROM ZERO
 122       IF(PIVOT) 3,15,3
 123 C        GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY
 124     3 ICOL=LPIV-IPIVC+1
 125       INDEX(MAIN)=ICOL
 126       IF(ICOL-MAIN) 6,6,4
 127 C        COMPLEMENT THE DETERMINANT
 128     4 DETER=-DETER
 129 C        POINTER TO LINE PIVOT FOUND
 130       ICOL=ICOL-IDIM
 131 C        POINTER TO EXACT PIVOT LINE
 132       I3=MAIN-IDIM
 133       DO 5 I=1,IEMAT
 134       ICOL=ICOL+IDIM
 135       I3=I3+IDIM
 136       SWAP=A(I3)
 137       A(I3)=A(ICOL)
 138     5 A(ICOL)=SWAP
 139 C        COMPUTE DETERMINANT
 140     6 DETER=DETER*PIVOT
 141       PIVOT=1./PIVOT
 142 C        TRANSFORM PIVOT COLUMN
 143       I3=IPIVC+NMIN1
 144       DO 7 I=IPIVC,I3
 145     7 A(I)=-A(I)*PIVOT
 146       A(IPIVC1)=PIVOT
 147 C        PIVOT ELEMENT TRANSFORMED
 148 C
 149 C        NOW CONVERT REST OF THE MATRIX
 150       I1=MAIN-IDIM
 151 C        POINTER TO PIVOT LINE ELEMENTS
 152       ICOL=1-IDIM
 153 C        GENERAL COLUMN POINTER
 154       DO 10 I=1,IEMAT
 155       ICOL=ICOL+IDIM
 156       I1=I1+IDIM
 157 C        POINTERS MOVED
 158       IF(I-MAIN) 8,10,8
 159 C        PIVOT COLUMN EXCLUDED
 160     8 JCOL=ICOL+NMIN1
 161       SWAP=A(I1)
 162       I3=IPIVC-1
 163       DO 9 I2=ICOL,JCOL
 164       I3=I3+1
 165     9 A(I2)=A(I2)+SWAP*A(I3)
 166       A(I1)=SWAP*PIVOT
 167    10 CONTINUE
 168    11 CONTINUE
 169 C        NOW REARRANGE THE MATRIX TO GET RIGHT INVERS
 170       DO 14 I1=1,N
 171       MAIN=N+1-I1
 172       LPIV=INDEX(MAIN)
 173       IF(LPIV-MAIN) 12,14,12
 174    12 ICOL=(LPIV-1)*IDIM+1
 175       JCOL=ICOL+NMIN1
 176       IPIVC=(MAIN-1)*IDIM+1-ICOL
 177       DO 13 I2=ICOL,JCOL
 178       I3=I2+IPIVC
 179       SWAP=A(I2)
 180       A(I2)=A(I3)
 181    13 A(I3)=SWAP
 182    14 CONTINUE
 183       DETERM=DETER
 184       NERROR=0
 185       RETURN
 186    15 NERROR=MAIN
 187       DETERM=DETER
 188       RETURN
 189       END
 190       SUBROUTINE D508R2 (A,DIM1,N1,DIM2,N2,INDEX,NERROR,DETERM)
 191       INTEGER DIM1,DIM,PIVCOL,PIVCO1,TOPX,ENDX,TOPCOL,ENDCOL,EMAT
 192       DIMENSION A(DIM1),INDEX(DIM1)
 193       DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP
 194       DIM=DIM1
 195       DETER=1.0D0
 196       N=N1
 197       EMAT=N+N2
 198       NMIN1=N-1
 199       PIVCOL=-DIM
 200 C     MAIN LOOP TO CREATE TRIANGULAR
 201       DO 10 MAIN=1,N
 202       PIVOT=0.0D0
 203       PIVCOL=PIVCOL+DIM+1
 204       PIVCO1=PIVCOL+N-MAIN
 205 C     SEARCH PIVOT
 206       DO 2 I1=PIVCOL,PIVCO1
 207       IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1
 208     1 PIVOT=A(I1)
 209       LPIV=I1
 210     2 CONTINUE
 211 C     IS PIVOT DIFFERENT FROM ZERO
 212       IF(PIVOT) 3,15,3
 213 C     IS IT NECESSARY TO BRING PIVOT TO DIAGONAL
 214     3 IF(LPIV-PIVCOL) 4,6,4
 215     4 DETER=-DETER
 216       LPIV=LPIV-DIM
 217       I1=PIVCOL-DIM
 218       DO 5 I2=MAIN,EMAT
 219       LPIV=LPIV+DIM
 220       I1=I1+DIM
 221       SWAP=A(I1)
 222       A(I1)=A(LPIV)
 223     5 A(LPIV)=SWAP
 224     6 DETER=DETER*PIVOT
 225       IF (MAIN .EQ. N)  GO TO 10
 226       PIVOT=1./PIVOT
 227 C     MODIFY PIVOT COLUMN
 228       I1=PIVCOL+1
 229       DO 7 I2=I1,PIVCO1
 230     7 A(I2)=A(I2)*PIVOT
 231 C     CONVERT THE SUBMATRIX AND RIGHT SIDES
 232       I3=PIVCOL
 233       IROW=MAIN+1
 234       DO 9 I1=IROW,N
 235       I3=I3+1
 236       I4=PIVCOL
 237       I5=I3
 238       DO 8 I2=IROW,EMAT
 239       I4=I4+DIM
 240       I5=I5+DIM
 241     8 A(I5)=A(I5)-A(I4)*A(I3)
 242     9 CONTINUE
 243    10 CONTINUE
 244       DETERM=DETER
 245       NERROR=0
 246 C     COMPUTE THE SOLUTIONS
 247       NO=N+1
 248       TOPX=NMIN1*DIM+1
 249       DO 13 I=NO,EMAT
 250       TOPX=TOPX+DIM
 251       ENDX=TOPX+N
 252       TOPCOL=N*DIM+1
 253       ENDCOL=TOPCOL+NMIN1
 254       DO 12 I1=1,NMIN1
 255       ENDX=ENDX-1
 256       TOPCOL=TOPCOL-DIM
 257       ENDCOL=ENDCOL-DIM-1
 258       A(ENDX)=A(ENDX)/A(ENDCOL+1)
 259       SWAP=A(ENDX)
 260       I3=TOPX-1
 261       DO 11 I2=TOPCOL,ENDCOL
 262       I3=I3+1
 263    11 A(I3)=A(I3)-A(I2)*SWAP
 264    12 CONTINUE
 265       A(TOPX)=A(TOPX)/A(1)
 266    13 CONTINUE
 267 C     LEFTADJUST THE SOLUTIONS
 268       I=-DIM
 269       TOPX=NMIN1*DIM+1
 270       ENDX=TOPX+NMIN1
 271       DO 14 I1=NO,EMAT
 272       TOPX=TOPX+DIM
 273       ENDX=ENDX+DIM
 274       I=I+DIM
 275       I3=I
 276       DO 14 I2=TOPX,ENDX
 277       I3=I3+1
 278    14 A(I3)=A(I2)
 279       RETURN
 280 C     ERROR EXIT
 281    15 NERROR=-1
 282       DETERM=DETER
 283       WRITE(6,100)MAIN,MAIN
 284       RETURN
 285   100 FORMAT(' LINEQ1 ..... THE ',I10,'. COLUMN OF THE MATRIX CONTAINS'
 286      1 ,' ONLY ZEROS AT THE',I10,'. ELIMINATIONSTEP')
 287       END