Logo Search packages:      
Sourcecode: design version File versions  Download package

mlmats.f

      FUNCTION isub(i,j)
C-----------------------------------------------------------------------------
C     Computes subscript in lower triangular matrix corresponding to (i,j)
C-----------------------------------------------------------------------------
      INTEGER i,j,isub
      IF(i-j)10,10,20
10    isub=i+j*(j-1)/2
      RETURN
20    isub=j+i*(i-1)/2
      RETURN
      END
        SUBROUTINE sqtria(vsq,vtri,n,k)
C----------------------------------------------------------------------------
C       k=1 : converts n x n square symmetric matrix vsq to lower triangular
C               form and stores result in vtri
C       k=2 : converts lower triangular matrix vtri to n x n uncompressed
C               square matrix
C       F. Harrell 6Sep90
C----------------------------------------------------------------------------
        DOUBLE PRECISION vsq(n,n),vtri(*)
        IF(k.EQ.1) THEN
                l=0
                DO i=1,n
                DO j=1,i
                l=l+1
                vtri(l)=vsq(i,j)
                END DO
                END DO
          ELSE
                DO i=1,n
                DO j=1,n
                vsq(i,j)=vtri(isub(i,j))
                END DO
                END DO
                ENDIF
        RETURN
        END
      subroutine inner(b,x,n,z)
C-----------------------------------------------------------------------------
C     Computes dot product of b and x, each of length n, returns result in z
C-----------------------------------------------------------------------------
      DOUBLE PRECISION b(1),x(1),z
      z=0D0
        DO i=1,n
        z=z+b(i)*x(i)
        end do
      return
      end
      SUBROUTINE SPROD(M,V,P,N)
C-----------------------------------------------------------------------------
C     MULTIPLIES N*N SYMMETRIC MATRIX M STORED IN COMPRESSED FORMAT BY
C     THE N*1 VECTOR V AND RETURNS THE N*1 VECTOR PRODUCT P
C-----------------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION M(1),V(1),P(1)
      DO 20 I=1,N
      PI=0D0
      II=I*(I-1)/2
      DO 10 J=1,N
      IF(I-J)2,4,4
    2 IR=I+J*(J-1)/2
      GO TO 10
    4 IR=J+II
   10 PI=PI+M(IR)*V(J)
      P(I)=PI
   20 CONTINUE
      RETURN
      END
      SUBROUTINE AVA(A,V,P,N)
C-----------------------------------------------------------------------------
C     V IS AN N X N SYMMETRIC MATRIX AND A IS AN N X 1 VECTOR.
C     THIS ROUTINE RETURNS P=A






















'VAC-----------------------------------------------------------------------------      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      DOUBLE PRECISION  A(1),V(1)      P=0D0      K=0      DO 10 I=1,N         AI=A(I)         DO 20 J=1,I            K=K+1            IF (I.EQ.J) THEN               P=P+AI*AI*V(K)            ELSE               P=P+2D0*AI*A(J)*V(K)            ENDIF20       CONTINUE10    CONTINUE      RETURN      END      SUBROUTINE avia(a,v,p,n,idx,nidx,nrank,eps,vsub,wv1,wv2,wv3,     &                  wv4,pivot)C----------------------------------------------------------------------------C       V is an n x n symmetric matrix and a is an n x 1 vector.C       Returns P=a' v**-1 a and nrank=rank(v), where 
C       a=a(idx(i),i=1,...,nidx), v=(v(idx(i),idx(i),i=1,...,nidx).
C       vsub is nidx x nidx scratch matrix and wv1-wv4 are scratch
C       vectors of length nidx (except for wv3 which is 2*nidx).  
C   pivot is scratch integer vector
C       of length nidx.  eps is singularity criterion, e.g. 1d-7.
C       Uses Fortran routines dqr (see S function qr) and dqrsl1 
C   (see S function solve).  In R these are dqrdc2 and dqrsl (args
C      differ too).
C
C       F. Harrell 20 Nov 90
C
C----------------------------------------------------------------------------
        DOUBLE PRECISION a(n),wv1(nidx),wv2(nidx),wv3(*),wv4(nidx),
     &          v(n,n),eps,vsub(nidx,nidx),p
        INTEGER idx(nidx),pivot(nidx),dim(2)
        k=nidx
C       CALL intpr("k",1,k,1)
        dim(1)=k
        dim(2)=k
                DO i=1,k
                wv4(i)=a(idx(i))
                pivot(i)=i
                        DO j=1,k
                        vsub(i,j)=v(idx(i),idx(j))
                        ENDDO
                ENDDO
C       CALL dblepr('wv4',3,wv4,k)
C       CALL dblepr('vsub',4,vsub,k*k)
        nrank=k
C        CALL dqr(vsub,dim,pivot,wv2,eps,wv3,nrank)
        CALL dqrdc2(vsub,dim,dim,dim,eps,nrank,wv2,pivot,wv3)
C       CALL intpr('nrank',5,nrank,1)
        IF(nrank.LT.k)RETURN
                DO i=1,k
                wv3(i)=wv4(i)
                ENDDO
        j=1
        i=100
C        CALL dqrsl1(vsub,dim,wv2,nrank,wv4,1,wv3,wv1,i,j)
        CALL dqrsl(vsub,dim,dim,nrank,wv2,wv4,wv3,wv1,wv1,
     &             wv3,wv3,i,j)
        p=0d0
                DO i=1,k
                p=p+wv4(i)*wv1(i)
                ENDDO
C       CALL intpr('dim',3,dim,2)
C       CALL dblepr('vsub',4,vsub,k*k)
C       CALL dblepr('wv1',3,wv1,k)
C       CALL dblepr('wv4',3,wv4,k)
C       CALL dblepr('p',1,p,1)
        RETURN
        END
      SUBROUTINE AVIA2(A,V,P,N,idx,nidx,nrank,eps,vsub,s,swept)
C----------------------------------------------------------------------------
C     V IS AN N X N SYMMETRIC square MATRIX AND A IS AN 
C     N X 1 VECTOR.
C     THIS ROUTINE RETURNS P=a







































































' vinverse a and nrank=rank(v) whereC     a=A(idx(i),i=1,...,nidx), v=V(idx(i),idx(i),i=1,...,nidx).C     S(nidx) is DOUBLE PRECISION scratch vector, SWEPT(nidx) is LOGICAL scratch C     vector, VSUB(nidx*(nidx+1)/2) is DOUBLE PRECISION scratch vectorC     eps is singularity criterion, e.g. 1D-6CC     F. Harrell 6 Sep90C----------------------------------------------------------------------------      IMPLICIT DOUBLE PRECISION (A-H,O-Z)      DOUBLE PRECISION  A(n),V(n,n),s(nidx),vsub(*)      INTEGER idx(nidx)      LOGICAL swept(nidx)        l=0                DO i=1,nidx                swept(i)=.FALSE.                idxi=idx(i)C       Initialize s vector to diagonal elements                s(i)=v(idxi,idxi)                DO j=1,i                l=l+1                vsub(l)=v(idxi,idx(j))                END DO                END DO        nrank=0                DO i=1,nidx                CALL GSWEEP(s,vsub,i,lsing,nidx,eps,swept,ifault)                IF(lsing.EQ.0)nrank=nrank+1                ENDDO      P=0D0      K=0      DO 10 I=1,NidxC       Singularities are like parameter never appeared        IF(swept(i)) THEN                 AI=A(idx(i))            ELSE                AI=0D0                ENDIF         DO 20 J=1,I            K=K+1            IF (I.EQ.J) THEN               P=P+AI*AI*Vsub(K)            ELSE               P=P+2D0*AI*A(idx(J))*Vsub(K)            ENDIF20       CONTINUE10    CONTINUEC       gsweep returns negative of inverse        P=-P      RETURN      END        SUBROUTINE ainvb(a, b, aib, k, tol, irank, pivot,      &                   wv1, wv2, wv3)C-----------------------------------------------------------------------C       Uses same Fortran subroutines as S function solve to accuratelyC       compute aib=a inverse * b, for k x k symmetric matrix a stored inC       lower triangular form and k x 1 vector b.  wv1(k,k), wv2(k), wv3(2*k)C       are DOUBLE PRECISION scratch arrays and pivot(k) is INTEGER scratch vector. C       tol is tolerance, e.g. 1d-7.C       IF irank (output) < k, result is not computed.  Index of singularC       column will be stored in pivot(k) if irank<k.C-----------------------------------------------------------------------        DOUBLE PRECISION a(1),b(k),aib(k),tol,wv1(k,k),wv2(k),wv3(*)        INTEGER pivot(k),dim(2)        CALL sqtria(wv1, a, k, 2)        dim(1)=k        dim(2)=k        DO i=1,k                pivot(i)=i                ENDDO        irank=kC        CALL dqr(wv1, dim, pivot, wv2, tol, wv3, irank)        CALL dqrdc2(wv1,dim,dim,dim,tol,irank,wv2,pivot,wv3)C       CALL dblepr('wv1
',3,wv1,10)C       CALL intpr('pivot
',5,pivot,k)C       CALL dblepr('wv2
',3,wv2,k)C       CALL dblepr('wv3
',3,wv3,k)C       CALL intpr('irank









',5,irank,1)        IF(irank.LT.k)RETURN        DO i=1,k                wv3(i)=b(i)                ENDDO        j=1        i=100C        CALL dqrsl1(wv1, dim, wv2, irank, b, 1, wv3, aib, i, j)        CALL dqrsl(wv1,dim,dim,irank,wv2,b,wv3,aib,aib,     &             wv3,wv3,i,j)C       CALL intpr('wv1
',3,wv1,10)C       CALL intpr('dim
',3,dim,2)C       CALL dblepr('wv2
',3,wv2,k)C       CALL intpr('irank
',5,irank,1)C       CALL dblepr('b
',1,b,k)C       CALL dblepr('wv3
',3,wv3,k)C       CALL dblepr('aib
',3,aib,k)C       CALL intpr('i
',1,i,1)C       CALL intpr('j












































































































































































































Generated by  Doxygen 1.6.0   Back to index