* Matrix subprograms * * The subprograms are in alphabetical order. * * GMPR print general matrix * SMAVAT product of general and symmetric matrix AVAT * SMINV symmetric matrix inversion and solution of equations * SMPR print symmetric matrix * SMPRV print (symmetric) covariance matrix and vector * SMTOS copy (part of) symmetric matrix to symmetric matrix SUBROUTINE GMPR(A,N,M,TEXT) * prints the N-by-M general matrix A * - - - ---- * CALL GMPR(A,N,M,TEXT) * REAL A(*) CHARACTER*(*) TEXT * * general matrix or vector * * ... IF(N.EQ.1.OR.M.EQ.1) THEN * special case of vector WRITE(*,101) N,M,TEXT WRITE(*,102) (A(L),L=1,N*M) ELSE * matrix WRITE(*,103) N,M,TEXT IJ=0 DO I=1,N WRITE(*,104) I,(A(IJ+L),L=1,M) IJ=IJ+M END DO END IF WRITE(*,102) 100 RETURN 101 FORMAT(1X,I4,' *',I3,' Vector ',A/) 102 FORMAT(10X,10G12.5) 103 FORMAT(1X,I4,' *',I3,' Matrix ',A/) 104 FORMAT(1X,I4,5X,10G12.5/(10X,10G12.5/)) END SUBROUTINE SMAVAT(V,A,W,N,M) * Multiply symmetric N-by-N matrix from the left with general M-by-N * matrix and from the right with the transposed of the same general * matrix to form symmetric M-by-M matrix (used for error * propagation). * - - - - * CALL SMAVA(V,A,W,N,M) * - * T * W = A * V * A * M*M M*N N*N N*M * * * where V = symmetric N-by-N matrix * A = general N-by-M matrix * W = symmetric M-by-M matrix * REAL V(1),A(1),W(1) DOUBLE PRECISION CIK * ... DO I=1,(M*M+M)/2 W(I)=0.0 END DO * IL=-N IJS=0 DO I=1,M IJS=IJS+I-1 IL=IL+N LKL=0 DO K=1,N CIK=0.0 LKL=LKL+K-1 LK=LKL DO L=1,K LK=LK+1 CIK=CIK+A(IL+L)*V(LK) END DO DO L=K+1,N LK=LK+L-1 CIK=CIK+A(IL+L)*V(LK) END DO JK=K IJ=IJS DO J=1,I IJ=IJ+1 W(IJ)=W(IJ)+CIK*A(JK) JK=JK+N END DO END DO END DO END SUBROUTINE SMINV(V,B,NARG,M,NRANK) * Obtain solution of a system of linear equations V * X = B with * symmetric matrix V and inverse (for M = 1) or matrix inversion * only (for M = 0). * * - - - - * CALL SMINV(V,B,N,M,NRANK) * - - ----- * * V = symmetric N-by-N matrix in symmetric storage mode * V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, . . . * replaced by inverse matrix * B = N-vector (for M = 0 use a dummy argument) * replaced by solution vector * M = see above * * * Method of solution is by elimination selecting the pivot on the * diagonal each stage. The rank of the matrix is returned in NRANK. * For NRANK ne N, all remaining rows and cols of the resulting * matrix V and the corresponding elements of B are set to zero. * SMINV can be used for a dimension up to 100 (see INVCDR). * REAL V(*),B(*),VKK,D,E COMMON/INVCDR/DR(2,100) DATA EPS/1.E-6/ * * construct table * N=IABS(NARG) IF(NARG.GT.0) THEN * reset flags DO I=1,N DR(1,I)=1.0 END DO NI=N ELSE * call with negative N - matrix partially inverted NI=0 DO I=1,N IF(DR(1,I).NE.0.0) NI=NI+1 END DO END IF II=0 DO I=1,N II=II+I DR(2,I)=ABS(V(II)) END DO * loop begin, with loop on all remaining rows/cols NRANK=N-NI DO I=1,NI * search for pivot and test for linearity and zero matrix K=0 JJ=0 KK=0 VKK=0.0 DO J=1,N JJ=JJ+J IF(DR(1,J).NE.0.0) THEN IF(ABS(V(JJ)).GT.VKK.AND.ABS(V(JJ)).GT.EPS*DR(2,J)) THEN VKK=ABS(V(JJ)) K=J KK=JJ END IF END IF END DO IF(K.EQ.0.OR.VKK.LT.1.0E-16) THEN * no pivot found, clear of matrix IJ=0 DO III=1,N IF(M.EQ.1.AND.DR(1,III).NE.0.0) B(III)=0.0 DO J=1,III IJ=IJ+1 IF(DR(1,III)+DR(1,J).NE.0.0) V(IJ)=0.0 V(IJ)=-V(IJ) END DO END DO GOTO 100 END IF * preparation for elimination NRANK=NRANK+1 DR(1,K)=0.0 D=1.0/V(KK) IF(ABS(D).LT.1.0E-16) D=0.0 V(KK)=-D IF(M.EQ.1) B(K)=B(K)*D c write(*,*) ' d and Bk ',D,B(K) JK=KK-K JL=0 * elimination DO J=1,N IF(J.EQ.K) THEN JK=KK JL=JL+J ELSE IF(J.LT.K) THEN JK=JK+1 ELSE JK=JK+J-1 END IF E =V(JK) IF(ABS(E).LT.1.0E-16) E=0.0 V(JK)=D*E IF(M.EQ.1) B(J) =B(J)-B(K)*E LK =KK-K DO L=1,J JL=JL+1 IF(L.EQ.K) THEN LK=KK ELSE IF(L.LT.K) THEN LK=LK+1 ELSE LK=LK+L-1 END IF V(JL)=V(JL)-V(LK)*E END IF END DO END IF END DO END DO * change sign IJ=0 DO I=1,N DO J=1,I IJ=IJ+1 V(IJ)=-V(IJ) END DO END DO 100 RETURN END SUBROUTINE SMPR(V,N,TEXT) * prints the symmetric n-by-n matrix v * ---- * CALL SMPRT(V,N,TEXT) * REAL V(*) CHARACTER*(*) TEXT * ... WRITE(*,101) N,N,TEXT IF(N.LE.0) GOTO 100 II=0 DO I=1,N WRITE(*,102) I,(V(II+J),J=1,I) II=II+I END DO WRITE(*,102) 100 RETURN 101 FORMAT(1X,I4,' by',I3,' symmetric matrix ',A/) 102 FORMAT(1X,I4,5X,10G12.5/(10X,10G12.5/)) END SUBROUTINE SMPRV(X,V,N) * Prints the n-vector X and the symmetric N-by-N covariance matrix * V, the latter as a correlation matrix. * * - - - * CALL DPRPAC(X,V,N) X = vector of parameters * V = covariance matrix * REAL X,V DIMENSION X(*),V(*) INTEGER MC(15) * WRITE(*,103) WRITE(*,101) II=0 DO I=1,N IJ=II II=II+I ERR=0.0 IF(V(II).GT.0.0) ERR=SQRT(V(II)) L=0 JJ=0 DO J=1,I JJ=JJ+J IJ=IJ+1 RHO=0.0 PD=V(II)*V(JJ) IF(PD.GT.0.0) RHO=V(IJ)/SQRT(PD) L=L+1 MC(L)=100.0*ABS(RHO)+0.5 IF(RHO.LT.0.0) MC(L)=-MC(L) IF(J.EQ.I.OR.L.EQ.15) THEN IF(J.LE.15) THEN IF(J.EQ.I) THEN WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L-1) ELSE WRITE(*,102) I,X(I),ERR,(MC(M),M=1,L) END IF ELSE IF(J.EQ.I) THEN WRITE(*,103) (MC(M),M=1,L-1) ELSE WRITE(*,103) (MC(M),M=1,L) END IF L=0 END IF END IF END DO END DO WRITE(*,104) 100 RETURN 101 FORMAT(9X,'Param',7X,'error',7X,'correlation coefficients'/) 102 FORMAT(1X,I5,2G12.4,1X,15I5) 103 FORMAT(31X,15I5) 104 FORMAT(33X,'(correlation coefficients in percent)') END SUBROUTINE SMTOS(V,I,W,J,N) * Copy symmetric N-by-N matrix or a N-by-N submatrix of a symmetric * matrix to another symmetric matrix * * - - - * CALL SMTOS(V,I,W,J,N) * - - * * N rows and columns of the matrix V, starting from diagonal element * (i,i), are copied to the matrix W, starting in diagonal element * (j,j). thus if a complete symmetric matrix has to be copied, i=j=1 * has to be used. * REAL V(*),W(*) * ... IM=(I*I+I)/2-1 JM=(J*J+J)/2-1 DO K=1,N DO L=1,K IM=IM+1 JM=JM+1 W(JM)=V(IM) END DO IM=IM+I-1 JM=JM+J-1 END DO END