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 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 00006240 C CHECK THAT SIZES ARE NOT GREATER THAN DIMENSIONS 00006250 C 00006260 IF(M.GT.IRDA.OR.N.GT.ICDA)GO TO 997 00006270 IF(N.GT.IRDB.OR.K.GT.ICDB)GO TO 997 00006280 IF(M.GT.IRDC.OR.K.GT.ICDC)GO TO 997 00006290 GO TO (10,40,80),ICODE 00006300 C 00006310 C CODING FOR A*B=C (ICODE=1) 00006320 C 00006330 10 DO 20 I=1,K 00006340 DO 20 J=1,M 00006350 20 C(J,I)=0.0D0 00006360 DO 30 I=1,K 00006370 DO 30 L=1,N 00006380 DO 30 J=1,M 00006390 30 C(J,I)=C(J,I)+A(J,L)*B(L,I) 00006400 C 00006420 C CODING FOR A*B=A (ICODE=2) 00006430 RETURN 00006410 C 00006440 40 IF(K.GT.ICDA)GO TO 998 00006450 DO 70 I=1,M 00006460 DO 60 L=1,K 00006470 SUM=0.0D0 00006480 DO 50 J=1,N 00006490 50 SUM=SUM+A(I,J)*B(J,L) 00006500 60 WORK(L)=SUM 00006510 DO 70 J=1,K 00006520 70 C(I,J)=WORK(J) 00006530 RETURN 00006540 C 00006550 C CODING FOR A*B=B (ICODE=3) 00006560 C 00006570 80 IF(M.GT.IRDB)GO TO 999 00006580 DO 110 I=1,K 00006590 DO 100 L=1,M 00006600 SUM=0.0D0 00006610 DO 90 J=1,N 00006620 90 SUM=SUM+A(L,J)*B(J,I) 00006630 100 WORK(L)=SUM 00006640 DO 110 J=1,M 00006650 110 C(J,I)=WORK(J) 00006660 RETURN 00006670 C 00006680 C WRITE ERROR MESSAGE AND STOP 00006690 C 00006700 997 WRITE (6,1000) STOP 00006720 998 WRITE(6,1010) STOP 00006740 999 WRITE(6,1020) STOP 00006760 1000 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF MATRIX > DIMENSION OF00006770 1 MATRIX') 00006780 1010 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF A*B > DIMENSION OF A 00006790 1(ICODE=2)') 00006800 1020 FORMAT('1',10X,'ERROR IN MATMUL ==> SIZE OF A*B > DIMENSION OF B 00006810 1(ICODE=3)') 00006820 END 00006830 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