SUBROUTINE COMBO1 (A,AT,IROWA,ICOLA,B,BT,IROWB,ICOLB,W,P,N,M,XCAP) C ********************************************************************* C * SUBROUTINE 'COMBO1' DOES A LEAST SQUARES (COMBINED * C * CASE) ADJUSTMENT USING ADDITIONAL CONSTRAINTS BETWEEN THE * C * UNKNOWN PARAMETERS TO OBTAIN UNIQUE VALUES FOR THE UNKNOWNS * C * PLUS THEIR VARIANCE-COVARIANCE MATRIX. COMBO1 CALCULATES THE * C * XCAP VECTOR. * C * INPUT: * C * A - IS THE DESIGN MATRIX FOR THE UNKNOWNS * C * AT - IS THE TRANSPOSE OF A * C * IROWA - IS THE ROW DIMENSION OF A * C * ICOLA - IS THE COLUMN DIMENSION OF A * C * B - IS THE DESIGN MATRIX FOR THE OBSERVABLES * C * BT - IS THE TRANSPOSE OF B * C * IROWB - IS THE ROW DIMENSION OF B * C * ICOLB - IS THE COLUMN DIMENSION OF B * C * W - IS THE VECTOR OF MISCLOSURE * C * P - IS THE WEIGHT MATRIX * C * N - IS AT * M INV * A ALL INVERSE * C * M - IS B * P INV * BT ALL INVERSE * C * XCAP - IS A VECTOR WHICH CONTAINS THE APPROXIMATIONS * C * FOR THE UNKNOWNS * C * * C * OUTPUT: * C * XCAP - THE SOLUTION VECTOR * C * N INV - THE INVERSE OF AT * M INV * A ALL INVERSE * C * * C * NOTE: THIS ROUTINE USES SYSTEM SUBROUTINES EXCEPT FOR CHOLD * C * AND MATMUL. * C * * C * BY: R. CALDWELL (DECEMBER, 1976) * C * * C ********************************************************************* C IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N ) REAL*8 M,N DIMENSION A(IROWA,ICOLA),AT(ICOLA,IROWA),B(IROWB,ICOLB),BT(ICOLB, @ IROWB),W(IROWA),M(IROWB,IROWB),N(ICOLA,ICOLA), @ XCAP(ICOLA),P(ICOLB,ICOLB) C C MULTIPLYING B * P INV BT C CALL TRNSD (BT,ICOLB,B,IROWB,IROWB,ICOLB) CALL TRNSD (AT,ICOLA,A,IROWA,IROWA,ICOLA) CALL MATMUL (B,IROWB,ICOLB,P,ICOLB,ICOLB,B,IROWB,ICOLB,IROWB, @ ICOLB,ICOLB,2) CALL MATMUL (B,IROWB,ICOLB,BT,ICOLB,IROWB,M,IROWB,IROWB,IROWB, @ ICOLB,IROWB,1) C C INVERTING M C CALL CHOLD (M,IROWB,IROWB,DETA,&100) C C MULTIPLYING AT * M INV * A C CALL MATMUL (AT,ICOLA,IROWA,M,IROWB,IROWB,AT,ICOLA,IROWA,ICOLA, @ IROWA,IROWA,2) CALL MATMUL (AT,ICOLA,IROWA,A,IROWA,ICOLA,N,ICOLA,ICOLA,ICOLA, @ IROWA,ICOLA,1) C C INVERTING N C CALL CHOLD (N,ICOLA,ICOLA,DETA,&100) C C MULTIPLYING AT * M INV * W C CALL MATMUL (AT,ICOLA,IROWA,W,IROWA,1,XCAP,ICOLA,1,ICOLA,IROWA,1, @ 1) C C DETERMINING THE SOLUTION VECTOR C CALL MATMUL (N,ICOLA,ICOLA,XCAP,ICOLA,1,XCAP,ICOLA,1,ICOLA,ICOLA, @ 1,3) C C RETURNING TO THE MAIN PROGRAM / SUBROUTINE C 100 RETURN END SUBROUTINE COMBO2 (A,IROWA,ICOLA,B,BT,IROWB,ICOLB,XCAP,M,N,W,P,DF, * V,VFACT,QX) C ********************************************************************* C * SUBROUTINE 'COMBO2' COMPLETES THE LEAST SQUARES * C * (COMBINED CASE) ADJUSTMENT AFTER THE SOLUTION VECTOR, XCAP HAS * C * BEEN SOLVED FOR. THIS SUBROUTINE COMPUTES THE RESIDUALS,THE * C * VARIANCE VACTOR AND THE VARIANCE - COVARIANCE MATRIX FOR THE * C * UNKNOWNS. * C * * C * INPUT: * C * A - IS THE DESIGN MATRIX FOR THE UNKNOWNS * C * IROWA - IS THE ROW DIMENSION OF A * C * ICOLA - IS THE COLUMN DIMENSION OF A * C * B - IS THE DESIGN MATRIX FOR THE OBSERVABLES * C * BT - IS THE TRANSPOSE OF B * C * IROWB - IS THE ROW DIMENSION OF B * C * ICOLB - IS THE COLUMN DIMENSION OF B * C * XCAP - IS A VECTOR WHICH CONTAINS THE APPROXIMATIONS * C * M - IS B * P INV * BT ALL INVERSE * C * N - IS AT * M INV * A ALL INVERSE * C * W - IS THE VECTOR OF MISCLOSURE * C * P - IS THE WEIGHT MATRIX * C * FOR THE UNKNOWNS * C * DF - THE DEGREES OF FREEDOM * C * * C * OUTPUT: * C * V - THE VECTOR OF RESIDUALS * C * VFACT - THE VARIANCE VACTOR * C * QX - THE VARIANCE - COVARIANCE MATRIX * C * * C * NOTE: THIS ROUTINE USES SYSTEM SUBROUTINES EXCEPT FOR CHOLD * C * AND MATMUL. * C * * C * BY: R. CALDWELL ( DECEMBER, 1979) * C * * C ********************************************************************* C IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N) REAL*8 M,N DIMENSION A(IROWA,ICOLA),B(IROWB,ICOLB),BT(ICOLB,IROWB), % XCAP(ICOLA),M(IROWB,IROWB),N(ICOLA,ICOLA),W(IROWA), % P(ICOLB,ICOLB),V(ICOLB),QX(ICOLA,ICOLA),VFACT(1,1) C C DETERMINATION OF THE CORRELATES C CALL MATMUL (A,IROWA,ICOLA,XCAP,ICOLA,1,V,ICOLB,1,IROWA,ICOLA,1,1) CALL MADDD (V,ICOLB,V,ICOLB,W,IROWA,IROWA,1) CALL MATMUL (M,IROWB,IROWB,V,ICOLB,1,V,ICOLB,1,IROWB,IROWB,1,3) C C DETERMINATION OF THE RESIDUALS C CALL TRNSD (BT,ICOLB,B,IROWB,IROWB,ICOLB) CALL MATMUL (P,ICOLB,ICOLB,BT,ICOLB,IROWB,BT,ICOLB,IROWB,ICOLB, @ ICOLB,IROWB,3) CALL MATMUL (BT,ICOLB,IROWB,V,ICOLB,1,V,ICOLB,1,ICOLB,IROWB,1,3) DO 10 I = 1,ICOLB 10 V(I) = -1.0D0 * V(I) C C DETERMINATION OF THE VARIANCE FACTOR C CALL TRNSD (BT,1,V,ICOLB,ICOLB,1) CALL MATMUL (BT,1,ICOLB,P,ICOLB,ICOLB,BT,1,ICOLB,1,ICOLB,ICOLB,2) CALL MATMUL (BT,1,ICOLB,V,ICOLB,1,VFACT,1,1,1,ICOLB,1,1) C VFACT(1,1) = VFACT(1,1) / DF C DO 20 I = 1,ICOLA DO 20 J = 1,ICOLA 20 QX(I,J) = N (I,J) * VFACT(1,1) RETURN END SUBROUTINE CHOLD (A,IRDA,NA,DETA,*) C ****************************************************************** C * * C * SUBROUTINE 'CHOLD' INVERTS A SQUARE MATRIX USING * C * THE CHOLESKY DECOMPOSITION ALGORITHUM. * C * * C * INPUT: * C * A - AN ARRAY CONTAINING A POSITIVE DEFINITE * C * SYMMETRIC INPUT MATRIX, * C * IRDA - THE ROW DIMENSION OF THE ARRAY * C * CONTAINING THE INPUT MATRIX, * C * NA- THE SIZE OF THE INPUT MATRIX, * C * OUTPUT: * C * A - CONTAINS THE INVERSE OF THE INPUT MATRIX * C * (INPUT DESTROYED) * C * DETA - IS THE DETERMINANT OF THE INPUT MATRIX * C * * - ERROR RETURN (TAKEN IF NA .LT. 1 OR IF * C * DETA .LT. SING) * C * * C * BY: D.E. WELLS (NOVEMBER, 1971) * C * MODIFIED BY: R. CALDWELL (OCTOBER, 1976)* C * * C ****************************************************************** C IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N) DIMENSION A(IRDA,NA) DATA SING/1D-40/ C C CHOLESKY DECOMPOSITION OF THE INPUT MATRIX INTO A TRIANGULAR MATRIX C IF( NA.LT.1 ) GOTO 180 DETA = A(1,1) A(1,1) = DSQRT(A(1,1)) IF( NA.EQ.1 ) GOTO 60 DO 10 I=2,NA A(I,1) = A(I,1) / A(1,1) 10 CONTINUE DO 50 J=2,NA SUM = 0.0D0 J1 = J - 1 DO 20 K = 1,J1 SUM = SUM + A(J,K) ** 2 20 CONTINUE DETA = DETA * (A(J,J) - SUM) A(J,J) = DSQRT( A(J,J) - SUM) IF( J.EQ.NA ) GOTO 50 J2 = J + 1 DO 40 I=J2,NA SUM = 0.0D0 DO 30 K=1,J1 SUM = SUM + A(I,K) * A(J,K) 30 CONTINUE A(I,J) = (A(I,J) - SUM) / A(J,J) 40 CONTINUE 50 CONTINUE 60 IF( DABS(DETA) .LT. SING ) GOTO 160 C C INVERSION OF LOWER TRIANGULAR MATRIX C DO 70 I=1,NA A(I,I) = 1.0D0 / A(I,I) 70 CONTINUE IF( NA .EQ. 1 ) GOTO 100 N1 = NA - 1 DO 90 J=1,N1 J2 = J + 1 DO 90 I = J2,NA SUM = 0.0D0 I1 = I - 1 DO 80 K = J,I1 SUM = SUM + A(I,K) * A(K,J) 80 CONTINUE A(I,J) = - A(I,I) * SUM 90 CONTINUE C C CONSTRUCTION OF THE INVERSE OF A C 100 DO 150 J = 1,NA IF( J.EQ.1 ) GOTO 120 J1 = J - 1 DO 110 I = 1,J1 A(I,J) = A(J,I) 110 CONTINUE 120 DO 140 I = J,NA SUM = 0.0D0 DO 130 K = I,NA SUM = SUM + A(K,J) * A(K,I) 130 CONTINUE A(I,J) = SUM 140 CONTINUE 150 CONTINUE RETURN 160 WRITE(6,170) DETA 170 FORMAT ('-',T3,'*** E R R O R *** SINGULAR MATRIX IN CHOLD. TH', * 'E DETERMINANT IS ',E20.5) RETURN 1 180 WRITE(6,190) 190 FORMAT ('-','*** E R R O R *** THE INPUT MATRIX IS OF DIMENSION', * ' ZERO IN CHOLD') RETURN 1 END SUBROUTINE MATMUL(A,IRDA,ICDA,B,IRDB,ICDB,C,IRDC,ICDC,M,N,K,ICODE) C ****************************************************************** C * * C * SUBROUTINE 'MATMUL' MULTIPLIES MATRICES COLUMN BY * C * COLUMN IN ONE OF THREE WAYS, A*B=C (ICODE=1) * C * A*B=A (ICODE=2) * C * A*B=B (ICODE=3) * C * FOR REAL*8 MATRICES UP TO 1000 * 1000 * C * * C * IN ALL CASES IT CHECKS THAT THE SIZE OF A MATRIX IS * C * NOT GREATER THAN ITS DIMENSION AND IF THE INDICATED * C * MULTIPLICATION IS POSSIBLE. * C * * C * * C * NOTE: IF ICODE .NE. 1 INPUT MATRICES ARE DESTROYED * C * * C * EXAMPLE: * C * ICODE = 1 CALL MATMUL (A,10,15,B,10,25,C,25,25,10,10,15,1) * C * ICODE = 2 CALL MATMUL (A,10,15,B,10,25,A,10,15,10,10,10,2) * C * ICODE = 3 CALL MATMUL (A,10,15,B,10,25,B,10,25,10,10,15,3) * C * * C * INPUT: * C * A,IRDA,ICDA - MATRIX A & ITS DECLARED DIMENSIONS * C * B,IRDB,ICDB - MATRIX B & ITS DECLARED DIMENSIONS * C * M*N - IS THE SIZE OF A * C * N*K - IS THE SIZE OF B * C * ICODE - THE TYPE OF MULTIPLICATION DESIRED SEE ABOVE* C * OUTPUT: * C * M*K - MATRIX C * C * * C * BY: C.A.M CHAMBERLAIN (FEBRUARY, 1975) * C * MODIFIED BY: R. CALDWELL (SEPTEMBER 1976)* C * * C ****************************************************************** C IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4 (I-N) DIMENSION A(IRDA,ICDA),B(IRDB,ICDB),C(IRDC,ICDC),WORK(1000) C C CHECK THAT THE SIZES ARE NOT GREATER THAN DIMENSIONS C IF( M .GT. IRDA .OR. N .GT. ICDA ) GOTO 120 IF( N .GT. IRDB .OR. K .GT. ICDB ) GOTO 120 IF( M .GT. IRDC .OR. K .GT. ICDC ) GOTO 120 GOTO (10,40,80),ICODE WRITE(6,1030) ICODE RETURN C C CODE FOR A*B=C ( ICODE=1 ) C 10 DO 20 I = 1,K DO 20 J = 1,M C(J,I) = 0.0D0 20 CONTINUE DO 30 I = 1,K DO 30 L = 1,N DO 30 J = 1,M C(J,I) = C(J,I) + A(J,L) * B(L,I) 30 CONTINUE RETURN C C CODING FOR A*B=A ( ICODE=2 ) C 40 IF( K .GT. ICDA ) GOTO 130 DO 70 I = 1,M DO 60 L = 1,K SUM = 0.0D0 DO 50 J = 1,N SUM = SUM + A(I,J) * B(J,L) 50 CONTINUE WORK(L) = SUM 60 CONTINUE DO 70 J = 1,K C(I,J) = WORK(J) 70 CONTINUE RETURN C C CODING FOR A*B=B ( ICODE=3 ) C 80 IF( M .GT. IRDB ) GOTO 140 DO 110 I = 1,K DO 100 L = 1,M SUM = 0.0D0 DO 90 J = 1,N SUM = SUM + A(L,J) * B(J,I) 90 CONTINUE WORK(L) = SUM 100 CONTINUE DO 110 J = 1,M C(J,I) = WORK(J) 110 CONTINUE RETURN C C WRITE ERROR MESSAGE AND RETURN C 120 WRITE(6,1000) RETURN 130 WRITE(6,1010) RETURN 140 WRITE(6,1020) RETURN 1000 FORMAT('-',T3,'*** E R R O R *** TERMINATION IN MATMUL DUE ', * 'TO THE SIZE OF AN INPUT MATRIX IS ',/,T25,'GREATER THAN', * ' THE ROW DIMENSION OF THE MATRIX') 1010 FORMAT ('-',T3,'*** E R R O R *** TERMINATION IN MATMUL DUE ', * 'THE SIZE OF A*B IS GREATER THAN THE DIMENSION OF A.') 1020 FORMAT ('-',T3,'*** E R R O R *** TERMINATION IN MATNUL DUE', * ' TO THE SIZE OF A*B IS GREATER THAN THE DIMENSION OF B.') 1030 FORMAT ('-',T3,'*** E R R O R *** TERMINATION IN MATMUL DUE ', * 'TO AN INVALID ICODE. ',I5,' WAS SPECIFIED BY THE USER.') END