* * $Id: linsq.F,v 1.1.1.1 1996/04/01 15:02:21 mclareni Exp $ * * $Log: linsq.F,v $ * Revision 1.1.1.1 1996/04/01 15:02:21 mclareni * Mathlib gen * * #include "gen/pilot.h" SUBROUTINE LINSQ(K,N,M,F,X,Y,W,DA,H,COV,QPRINT,QMIN) DIMENSION X(M,N),F(K),Y(N),W(N),H(K,K),DA(K) DIMENSION INDEX(100) DOUBLE PRECISION H,DETERM DOUBLE PRECISION DWI,DYI,DFK1 KPLUS=K+1 APRINT=QPRINT DO 11 I=1,K DO 11 J=1,KPLUS 11 H(I,J) = 0.D0 IFLAG=1 DO 2 I=1,N CALL FCN(K,M,F,X(1,I),IFLAG) IFLAG=4 DWI = DBLE(W(I)) DYI = DBLE(Y(I)) DO 2 K1=1,K DFK1 = DBLE(F(K1)) H(K1,K+1) = H(K1,K+1) + DWI*DFK1*DYI DO 2 L=1,K1 2 H(K1,L) = H(K1,L) + DWI*DFK1*DBLE(F(L)) DO 3 J=1,K DO 3 I=1,J 3 H(I,J)=H(J,I) IF(COV)5,4,5 4 CALL D508R2(H,K,K,K+1,1,INDEX,NERROR,DETERM) MM=1 QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM) DO 7 I=1,K 7 F(I)=H(I,1) IF(APRINT.NE.0.) GO TO 10 RETURN 5 CALL D508R1(H,K,K,K+1,1,INDEX,NERROR,DETERM) MM=K+1 QMIN = QLINSQ(K,N,M,F,X,Y,W,DA,H,MM) DO 6 I=1,K 6 F(I)=H(I,K+1) DO 8 I=1,K HDA=H(I,I) 8 DA(I)=SQRT(HDA) IF(APRINT.EQ.0.) RETURN C C PRINTS ERRORS IN THE COEFFICIENTS C WRITE(6,100) 100 FORMAT(45X,'ERRORS IN THE COEFFICIENTS') WRITE(6,101)(DA(I),I=1,K) 101 FORMAT(5X,5G16.6) C C PRINTS LOWER DIAGONAL COVARIANCE MATRIX C WRITE(6,102) 102 FORMAT(///40X,'VARIANCE-COVARIANCE MATRIX'///) DO 9 I=1,K DO 12 J=1,K 12 F(J)=H(I,J) 9 WRITE(6,103)(F(JJ),JJ=1,I) 103 FORMAT(10X,5G20.6//) DO 13 I=1,K 13 F(I)=H(I,K+1) 10 CONTINUE C C PRINTS COEFFICIENTS AND SUM OF SQUARES C WRITE(6,104)QMIN 104 FORMAT(///30X,'SUM OF SQUARES=',G20.6///) WRITE(6,105) 105 FORMAT(45X,'COEFFICIENTS'///) WRITE(6,103)(F(I),I=1,K) RETURN END FUNCTION QLINSQ(K,N,M,F,X,Y,W,DA,H,MM) C COMPUTES THE SUM OF SQUARES DIMENSION F(K),X(M,N),DA(K),H(K,K),Y(N),W(N) DOUBLE PRECISION H IFLAG = 4 Q=0. DO 3 I=1,N CALL FCN(K,M,F,X(1,I),IFLAG) HSUM=0. DO 2 I1=1,K 2 HSUM=H(I1,MM)*F(I1)+HSUM 3 Q=(HSUM-Y(I))*(HSUM-Y(I))*W(I)+Q QLINSQ = Q RETURN END SUBROUTINE D508R1 (A,IDIM1,N1,IDIM2,N2,INDEX,NERROR,DETERM) C C MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP DIMENSION A(IDIM1),INDEX(IDIM1) DETER=1.0D0 N=N1 IEMAT=N+N2 IDIM=IDIM1 NMIN1=N-1 C THE ROUTINE DOES ITS OWN EVALUATION FOR DOUBLE SUBSCRIPTING OF C ARRAY A. IPIVC=1-IDIM C MAIN LOOP TO INVERT THE MATRIX DO 11 MAIN=1,N PIVOT=0.0D0 IPIVC=IPIVC+IDIM C SEARCH FOR NEXT PIVOT IN COLUMN MAIN. IPIVC1=IPIVC+MAIN-1 IPIVC2=IPIVC +NMIN1 DO 2 I1=IPIVC1,IPIVC2 IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1 1 PIVOT=A(I1) LPIV=I1 2 CONTINUE C IS PIVOT DIFFERENT FROM ZERO IF(PIVOT) 3,15,3 C GET THE PIVOT-LINE INDICATOR AND SWAP LINES IF NECESSARY 3 ICOL=LPIV-IPIVC+1 INDEX(MAIN)=ICOL IF(ICOL-MAIN) 6,6,4 C COMPLEMENT THE DETERMINANT 4 DETER=-DETER C POINTER TO LINE PIVOT FOUND ICOL=ICOL-IDIM C POINTER TO EXACT PIVOT LINE I3=MAIN-IDIM DO 5 I=1,IEMAT ICOL=ICOL+IDIM I3=I3+IDIM SWAP=A(I3) A(I3)=A(ICOL) 5 A(ICOL)=SWAP C COMPUTE DETERMINANT 6 DETER=DETER*PIVOT PIVOT=1./PIVOT C TRANSFORM PIVOT COLUMN I3=IPIVC+NMIN1 DO 7 I=IPIVC,I3 7 A(I)=-A(I)*PIVOT A(IPIVC1)=PIVOT C PIVOT ELEMENT TRANSFORMED C C NOW CONVERT REST OF THE MATRIX I1=MAIN-IDIM C POINTER TO PIVOT LINE ELEMENTS ICOL=1-IDIM C GENERAL COLUMN POINTER DO 10 I=1,IEMAT ICOL=ICOL+IDIM I1=I1+IDIM C POINTERS MOVED IF(I-MAIN) 8,10,8 C PIVOT COLUMN EXCLUDED 8 JCOL=ICOL+NMIN1 SWAP=A(I1) I3=IPIVC-1 DO 9 I2=ICOL,JCOL I3=I3+1 9 A(I2)=A(I2)+SWAP*A(I3) A(I1)=SWAP*PIVOT 10 CONTINUE 11 CONTINUE C NOW REARRANGE THE MATRIX TO GET RIGHT INVERS DO 14 I1=1,N MAIN=N+1-I1 LPIV=INDEX(MAIN) IF(LPIV-MAIN) 12,14,12 12 ICOL=(LPIV-1)*IDIM+1 JCOL=ICOL+NMIN1 IPIVC=(MAIN-1)*IDIM+1-ICOL DO 13 I2=ICOL,JCOL I3=I2+IPIVC SWAP=A(I2) A(I2)=A(I3) 13 A(I3)=SWAP 14 CONTINUE DETERM=DETER NERROR=0 RETURN 15 NERROR=MAIN DETERM=DETER RETURN END SUBROUTINE D508R2 (A,DIM1,N1,DIM2,N2,INDEX,NERROR,DETERM) INTEGER DIM1,DIM,PIVCOL,PIVCO1,TOPX,ENDX,TOPCOL,ENDCOL,EMAT DIMENSION A(DIM1),INDEX(DIM1) DOUBLE PRECISION A,DETERM,DETER,PIVOT,SWAP DIM=DIM1 DETER=1.0D0 N=N1 EMAT=N+N2 NMIN1=N-1 PIVCOL=-DIM C MAIN LOOP TO CREATE TRIANGULAR DO 10 MAIN=1,N PIVOT=0.0D0 PIVCOL=PIVCOL+DIM+1 PIVCO1=PIVCOL+N-MAIN C SEARCH PIVOT DO 2 I1=PIVCOL,PIVCO1 IF(ABS(A(I1))-ABS(PIVOT)) 2,2,1 1 PIVOT=A(I1) LPIV=I1 2 CONTINUE C IS PIVOT DIFFERENT FROM ZERO IF(PIVOT) 3,15,3 C IS IT NECESSARY TO BRING PIVOT TO DIAGONAL 3 IF(LPIV-PIVCOL) 4,6,4 4 DETER=-DETER LPIV=LPIV-DIM I1=PIVCOL-DIM DO 5 I2=MAIN,EMAT LPIV=LPIV+DIM I1=I1+DIM SWAP=A(I1) A(I1)=A(LPIV) 5 A(LPIV)=SWAP 6 DETER=DETER*PIVOT IF (MAIN .EQ. N) GO TO 10 PIVOT=1./PIVOT C MODIFY PIVOT COLUMN I1=PIVCOL+1 DO 7 I2=I1,PIVCO1 7 A(I2)=A(I2)*PIVOT C CONVERT THE SUBMATRIX AND RIGHT SIDES I3=PIVCOL IROW=MAIN+1 DO 9 I1=IROW,N I3=I3+1 I4=PIVCOL I5=I3 DO 8 I2=IROW,EMAT I4=I4+DIM I5=I5+DIM 8 A(I5)=A(I5)-A(I4)*A(I3) 9 CONTINUE 10 CONTINUE DETERM=DETER NERROR=0 C COMPUTE THE SOLUTIONS NO=N+1 TOPX=NMIN1*DIM+1 DO 13 I=NO,EMAT TOPX=TOPX+DIM ENDX=TOPX+N TOPCOL=N*DIM+1 ENDCOL=TOPCOL+NMIN1 DO 12 I1=1,NMIN1 ENDX=ENDX-1 TOPCOL=TOPCOL-DIM ENDCOL=ENDCOL-DIM-1 A(ENDX)=A(ENDX)/A(ENDCOL+1) SWAP=A(ENDX) I3=TOPX-1 DO 11 I2=TOPCOL,ENDCOL I3=I3+1 11 A(I3)=A(I3)-A(I2)*SWAP 12 CONTINUE A(TOPX)=A(TOPX)/A(1) 13 CONTINUE C LEFTADJUST THE SOLUTIONS I=-DIM TOPX=NMIN1*DIM+1 ENDX=TOPX+NMIN1 DO 14 I1=NO,EMAT TOPX=TOPX+DIM ENDX=ENDX+DIM I=I+DIM I3=I DO 14 I2=TOPX,ENDX I3=I3+1 14 A(I3)=A(I2) RETURN C ERROR EXIT 15 NERROR=-1 DETERM=DETER WRITE(6,100)MAIN,MAIN RETURN 100 FORMAT(' LINEQ1 ..... THE ',I10,'. COLUMN OF THE MATRIX CONTAINS' 1 ,' ONLY ZEROS AT THE',I10,'. ELIMINATIONSTEP') END