C #LSFTB V1B 8-JUN-73 C LAST UPDATE C SUBROUTINE LSFTB(DT,WT,FT,M,ST,ER,CO,N,N1,CHISQ) DIMENSION DT(M), WT(M), FT(M), ST(M,N), ER(N1,N), CO(N) C C *PURPOSE. C TO PERFORM A GENERAL LINEAR LEAST SQUARES FIT. C C *PARAMETERS. C DT - IS A REAL ARRAY OF DATA TO BE FITTED OF LENGTH "M" C WT - IS A REAL ARRAY OF WEIGHTS FOR THE DATA OF LENGTH "M" C FT - IS A REAL ARRAY RETURNED WITH VALUES OF THE FIT, LENGTH "M" C M - IS THE NUMBER OF DATA POINTS. C ST - A REAL MATRIX CONTAINING THE COMPONENTS TO BE FITTED TO THE C DATA. THE COMPONENTS ARE STORED IN THE COLUMNS OF "ST" C AND THE FIRST DIMENSION OF "ST" MUST BE "M". THE NUMBER C OF COMPONENTS FITTED IS GIVEN BY "N", AND SO THE SECOND C DIMENSION OF "ST" MUST BE AT LEAST "N". C ER - IS RETURNED AS THE INVERSE OF THE NORMAL MATRIX. I.E. ITS C DIAGONAL ELEMENTS ARE THE THEORETICAL VARIANCES OF THE C COEFFICIENTS OF THE FIT "CO" WHEN MULTIPLIED BY "CHISQ". C ITS OFF-DIAGONAL ELEMENTS ARE THE COVARIANCES BETWEEN PAIRS C OF COEFFICIENTS. THE FIRST DIMENSION OF "ER" MUST BE N1 C AND THE SECOND AT LEAST N. C CO - IS A REAL ARRAY OF LENGTH "N" OR MORE WHICH IS RETURNED WITH C THE COEFFICIENTS OF THE FIT. C N - IS THE NUMBER OF COMPONENTS BEING FITTED. C N1 - IS THE FIRST DIMENSION OF "ER". I.E. THE MAX NUMBER OF C COMPONENTS ALLOWED FOR IN "ER". C CHISQ - IS THE "NUCLEAR PHYSICIST'S CHI-SQUARE!" . THAT IS THE SUM C OF THE SQUARES OF THE DIFFERENCES BETWEEN FIT AND DATA, C MULTIPLIED BY THE WEIGHT FOR THE DATA, ALL DIVIDED BY (M-N) C THE NUMBER OF DEGREES OF FREEDOM. FOR A REASONABLE FIT AND C SENSIBLE WEIGHTS(ERRORS) THIS VALUE SHOULD BE ABOUT UNITY. C C INPUT VARIABLES: DT, WT, M, ST, N, N1, CHISQ C OUTPUT VARIABLES: FT, ER, CO, CHISQ C C *METHOD. C THE NORMAL MATRIX IS BUILT UP IN THE STRAIGHTFORWARD WAY, IT IS C INVERTED BY CALLING "INVTB" FROM THE LIBRARY THIS BEING A C GAUSSIAN ELIMINATION TYPE OF INVERSION. THE CALCULATION OF COEFFICIENTS C AND FIT AND CHISQ IS STRAIGHTFORWARD. THE CODE HAS BEEN PARTIALLY C OPTIMIZED TO REDUCE THE NUMBER OF FLOATING POINT MULTIPLICATIONS C TO A MINIMUM AT THE EXPENSE OF ARRAY INDEXING AND CODE LENGTH. C NOTE THERE IS A SPECIAL ENTRY FEATURE. AFTER CALLING C THE ROUTINE ONCE, A SUBSEQUENT CALL IN WHICH THE ONLY INPUT C PARAMETER CHANGED IS "DT" CAN BE SPEEDED UP CONSIDERABLY BY SETTING C CHISQ -VE. OTHERWISE CHISQ SHOULD NOT BE SET -VE! C C *ACCURACY. C THE ACCURACY CAN BE AFFECTED BY THE CHARACTER OF THE COMPONENTS C BEING FITTED. THE IDEAL SITUATION IS THAT ALL COMPONENTS SHOULD C BE ORTHOGONAL AND THE WORST SITUATION IS WHERE SOME ARE PARALLEL C IN WHICH CASE THEIR INDIVIDUAL COEFFICIENTS CAN NOT BE DETERMINED C PROPERLY. IN PRACTICE THE USER SHOULD TRY TO AVOID COMPONENTS WHICH C ARE VERY SIMILAR IN SHAPE OR WHICH PRODUCE LARGE OFF-DIAGONAL C ELEMENTS IN "ER". SENSIBLY USED THE ROUTINE SHOULD GIVE AT LEAST C FOUR SIGNIFICANT DIGITS. NOTE THIS REFERS TO MACHINE ROUNDING AND C TRUNCATION ERROR ONLY!! C C *RESTRICTIONS. C "LSFTB" DOES NOT CHECK IN ANY WAY THAT THE NORMAL MATRIX C CAN BE INVERTED SUCCESSFULLY. THE CALLER SHOULD THEREFORE BE C CAREFUL NOT TO USE COMPONENTS OF IDENTICAL OR NEARLY SO SHAPES. C C *ERROR CONDITIONS. C C *NON STANDARD ROUTINES CALLED. C INVTB (MRMLIB). C C *TYPICAL TIMES. C C *ORIGIN. M.R.MANNING. C C *COMMENTS. C "LSFTB" HAS BEEN DERIVED FROM "LSFTA". THE ONLY DIFFERENCE C BEING THAT "LSFTB" USES A MATRIX INVERSION ROUTINE WITH C COLUMN PIVOTING. WHILE PIVOTING IS NOT STRICTLY NECESSARY C TO INVERT THE POSITIVE DEFINITE MATRICES INVOLVED HERE IT C HAS BEEN FOUND THAT PIVOTING DOES REDUCE THE ROUNDING C ERROR SIGNIFICANTLY PARTICULARLY WHEN THE COMPONENTS C BEING FITTED ARE SIMILAR. C C #END. C C #START. C C CHECK ENTRY TYPE: IF(CHISQ) 50,10,10 C C NORMAL ENTRY: COMPOSE AND INVERT NORMAL MATRIX C C CLEAR MATRIX 10 DO 15 K=1,N DO 15 J=1,K 15 ER(K,J)=0. C BUILD UP MATRIX DATA POINT BY DATA POINT. DO 20 L=1,M WL=WT(L) DO 20 K=1,N XK=ST(L,K)*WL DO 20 J=1,K 20 ER(K,J)=ER(K,J)+XK*ST(L,J) C C THAT FIXES LOWER TRIANGLE: NOW SYMMETRIZE. C DO 30 K=1,N DO 30 J=1,K 30 ER(J,K)=ER(K,J) C C INVERT NORMAL MATRIX IN SITU WITH STANDARD ROUTINE. CALL INVTB(ER,N,N1,DET, CO) C NOTE CO IS USED AS A WORKING ARRAY BY "INVTB" C C PHASE 2: SPECIAL ENTRY: COMPOSE RIGHT HAND SIDE OF C NORMAL EQUATIONS AND MULTIPLY THROUGH BY INVERSE NORMAL MATRIX. C C CLEAR COEFFICIENT ARRAY 50 DO 60 J=1,N 60 CO(J)=0. C C CALCULATE RIGHT HAND SIDE SIGMA[ ST(L,K)*WT(L)*DT(L)] C AND MULTIPLY BY INVERTED MATRIX ER(J,K) ALL AT ONCE. C DO IT DATA POINT BY DATA POINT. C DO 70 L=1,M WL=WT(L)*DT(L) DO 70 K=1,N XK=ST(L,K)*WL DO 70 J=1,N 70 CO(J)=CO(J) + XK*ER(J,K) C C NOW FINISH OFF BY CALCULATING FIT AT EACH POINT AND THE C TOTAL CHISQUARED. C CHISQ=0. DO 80 L=1,M XK=0. WL=WT(L) DO 75 K=1,N 75 XK=XK+ST(L,K)*CO(K) FT(L)=XK XK=XK-DT(L) 80 CHISQ=CHISQ+WL*XK*XK C ENSURE DEGREES OF FREEDOM NOT ZERO. 95 CHISQ=CHISQ/AMAX0(M-N, 1) C END