SUBROUTINE LSA (LU,F,PF,A,PX, $ NF,IRPF,ICPF,NX,IRPX,ICPX,IRA,IRAPA, $ X,APA,R,RNORM,APB,DET,IERR) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION A(IRA,NX),APA(IRAPA,NX),APB(NX), $ F(NF),PF(IRPF,ICPF),PX(IRPX,ICPX), $ R(NF),X(NX) C$ C NAME LSA C C TYPE SUBROUTINE C C PURPOSE LSA = LEAST SQUARES APPROXIMATION C SOLVES EQUATION X = (PX + A'.PF.A)" .A'.PF.F C WHERE ' MEANS TRANSPOSE C " MEANS INVERSE C C AUTHOR/DATE D.WELLS 780315 C C EXTERNALS SPIN C C IRAPA,X,APA,R,RNORM,APB,DET,IERR) C C CALLING CALL LSA(LU,F,PF,A,PX,NF,IRPF,ICPF,NX,IRPX,ICPX,IRA, C INPUT PARAMETERS C LU = LISTING LU C F(NF) = FUNCTION TO BE APPROXIMATED (OBSERVATIONS) C PF(NF,NF) = WEIGHT MATRIX OF F C A(NF,NX) = DESIGN MATRIX C PX(NX,NX) = APRIORI WEIGHT MATRIX OF X C NF = NUMBER OF OBSERVATIONS C IRPF = PF ROW DIMENSION IN CALLING PROGRAM C ICPF = PF COLUMN DIMENSION C NX = NUMBER OF UNKNOWNS C IRPX = PX ROW DIMENSION IN CALLING PROGRAM C ICPX = PX COLUMN DIMENSION C IRA = A ROW DIMENSION IN CALLING ROUTINE C IRAPA = APA ROW DIMENSION IN CALLING ROUTINE C C OUTPUT PARAMETERS C X(NX) = APPROXIMATING COEFFICIENTS (UNKNOWNS) C APA(NX,NX) = (PX + A'.PF.A')" = RELATIVE COVARIANCE OF X C R(NF) = P(NF) - F(NF) = RESIDUALS C RNORM = QUADRATIC NORM OF R C APB(NX) = A'.PF.F = NORMAL EQUATION VECTOR C DET = DETERMINANT OF APA" C IERR = 0 SUCCESSFUL RETURN C 1 (PX + A'.PF.A) IS SINGULAR IN SPIN C 2 NX IS ZERO IN SPIN C 3 NF .LT. NX AND PX IS NULL C 4 IRPF MUST BE 1 OR .GE. NF C 5 ICPF MUST BE 1 OR NF C 6 IRPX MUST BE 1 OR .GE. NX C 7 ICPX MUST BE 1 OR NX C 8 IRA MUST BE .GE. NF C 9 IRAPA MUST BE .GE. NX C C LANGUAGE FORTRAN C C COMMENTS INPUT PF CAN BE IDENTITY (IRPF = 1, ICPF = 1) C DIAGONAL (IRPF = 1, ICPF = NF) C FULL (IRPF .GE. NF, ICPF = NF) C INPUT PX CAN BE NULL (IRPX = 1, ICPX = 1) C DIAGONAL (IRPX = 1, ICPX = NX) C FULL (IRPX .GE. NX, ICPX = NX) C$ C CHECK INPUT DIMENSIONS IERR = 0 IF(NF .LT. NX .AND. ICPX .NE. NX) IERR = 3 IF(IRPF .NE. 1 .AND. IRPF .LT. NF) IERR = 4 IF(ICPF .NE. 1 .AND. ICPF .NE. NF) IERR = 5 IF(IRPX .NE. 1 .AND. IRPX .LT. NX) IERR = 6 IF(ICPX .NE. 1 .AND. ICPX .NE. NX) IERR = 7 IF(IRA .LT. NF) IERR = 8 IF(IRAPA .LT. NX) IERR = 9 IF(IERR .NE. 0) RETURN C CLEAR ARRAYS APA AND APB DO 5 I =1,NX APB(I) = 0. DO 5 J = 1,NX 5 APA(I,J) = 0. RNORM = 0. C COMPUTE APA = A'.PF.A AND APB = A'.PF.F IF(ICPF .EQ. NF) GO TO 15 C PF = IDENTITY MATRIX DO 10 K=1,NF DO 10 J = 1,NX APB(J) = APB(J) + A(K,J) * F(K) DO 10 I = 1,NX 10 APA(I,J) = APA(I,J) + A(K,J) * A(K,I) GO TO 35 15 IF(IRPF .NE. 1) GO TO 25 C PF = DIAGONAL MATRIX DO 20 K = 1,NF DO 20 J = 1,NX APB(J) = APB(J) + A(K,J) * PF(1,K) * F(K) DO 20 I = 1,NX 20 APA(I,J) = APA(I,J) + A(K,J) * PF(1,K) * A(K,I) GO TO 35 C PF = FULL MATRIX 25 DO 30 K = 1,NF DO 30 J = 1,NX DO 30 L = 1,NF APB(J) = APB(J) + A(K,J) * PF(K,L) * F(L) DO 30 I = 1,NX 30 APA(I,J) = APA(I,J) + A(K,J) * PF(K,L) * A(L,I) C ADD PX TO APA 35 IF(ICPX .NE. NX) GO TO 55 IF(IRPX .NE. 1) GO TO 45 C PX = DIAGONAL MATRIX DO 40 I = 1,NX 40 APA(I,I) = APA(I,I) + PX(1,I) GO TO 55 C PX = FULL MATRIX 45 DO 50 I = 1,NX DO 50 J = I,NX 50 APA(J,I) = APA(J,I) + PX(J,I) C CHOLESKI INVERSION OF APA 55 CALL SPIN(APA,NX,IRAPA,DET,IDEXP) IF(IERR .NE. 0) RETURN C COMPUTE COEFFICIENTS X = APA * APB DO 60 I = 1,NX X(I) = 0. DO 60 J = 1,NX 60 X(I)= X(I) + APA(I,J) * APB(J) C COMPUTE APPROXIMANT (P = A * X) AND RESIDUALS (R = P - F) DO 70 I = 1,NF P = 0. DO 65 J = 1,NX 65 P = P + A(I,J) * X(J) 70 R(I)= P - F(I) C COMPUTE QUADRATIC NORM OF R IF(ICPF .EQ. NF) GO TO 80 C PF = IDENTITY MATRIX DO 75 I = 1,NF 75 RNORM = RNORM + R(I)**2 RETURN 80 IF(IRPF .NE. 1) GO TO 90 C PF = DIAGONAL MATRIX DO 85 I = 1,NF 85 RNORM = RNORM + R(I) ** 2 * PF(1,I) RETURN C PF = FULL MATRIX 90 DO 95 I = 1,NF DO 95 J = 1,NF 95 RNORM = RNORM + R(I) * PF(I,J) * R(J) RETURN END