+c-----------------------------------------------------------------------------
+ SUBROUTINE MATOUTR(N,A)
+c Prints the lower fragment of a symmetric matix
+ implicit none
+ integer n
+ double precision a(n*(n+1)/2)
+ integer i,j,k,nlim,jlim,jlim1
+ CHARACTER*6 LINE6 / '------' /
+ CHARACTER*12 LINE12 / '------------' /
+ double precision B(10)
+ include 'COMMON.IOUNITS'
+ DO 1 I=1,N,10
+ NLIM=MIN0(I+9,N)
+ WRITE (IOUT,1000) (K,K=I,NLIM)
+ WRITE (IOUT,1020) LINE6,(LINE12,K=I,NLIM)
+ 1000 FORMAT (/7X,10(I5,2X))
+ 1020 FORMAT (A6,10A7)
+ DO 2 J=I,N
+ JLIM=MIN0(J,NLIM)
+ JLIM1=JLIM-I+1
+ DO 3 K=I,JLIM
+ 3 B(K-I+1)=A(J*(J-1)/2+K)
+ WRITE (IOUT,1010) J,(B(K),K=1,JLIM1)
+ 2 CONTINUE
+ 1 CONTINUE
+ 1010 FORMAT (I3,3X,10(F7.2))
+ RETURN
+ END
+#ifdef FIVEDIAG
+c---------------------------------------------------------------------------
+ subroutine fivediagmult(n,DM,DU1,DU2,x,y)
+ implicit none
+ integer n
+ double precision DM(n),DU1(n),DU2(n),x(n),y(n)
+ integer i
+ y(1)=DM(1)*x(1)+DU1(1)*x(2)+DU2(1)*x(3)
+ y(2)=DU1(1)*x(1)+DM(2)*x(2)+DU1(2)*x(3)+DU2(2)*x(4)
+ do i=3,n-2
+ y(i)=DU2(i-2)*x(i-2)+DU1(i-1)*x(i-1)+DM(i)*x(i)
+ & +DU1(i)*x(i+1)+DU2(i)*x(i+2)
+ enddo
+ y(n-1)=DU2(n-3)*x(n-3)+DU1(n-2)*x(n-2)+DM(n-1)*x(n-1)
+ & +DU1(n-1)*x(n)
+ y(n)=DU2(n-2)*x(n-2)+DU1(n-1)*x(n-1)+DM(n)*x(n)
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine fivediaginv_mult(ndim,forces,d_a_vec)
+ implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.LAGRANGE.5diag'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ integer ndim
+ double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
+ & xsolv(ndim),d_a_vec(6*nres)
+ integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
+#ifdef TIMING
+ include 'COMMON.TIME1'
+ double precision time01
+ time01=MPI_Wtime()
+#endif
+ accel(:,:2*nres)=0.0d0
+ do j=1,3
+Compute accelerations in Calpha and SC
+ do ichain=1,nchain
+ n=dimen_chain(ichain)
+ iposc=iposd_chain(ichain)
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ do i=iposc,iposc+n-1
+ rs(i-iposc+1)=forces(3*(i-1)+j)
+ enddo
+#ifdef DEBUG
+ write (iout,*) "j",j," chain",ichain
+ write (iout,*) "rs"
+ write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
+ call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
+#ifdef DEBUG
+ write (iout,*) "xsolv"
+ write (iout,'(f10.5)') (xsolv(i),i=1,n)
+#endif
+ ind=1
+ do i=innt,inct
+ if (itype(i).eq.10)then
+ accel(j,i)=xsolv(ind)
+ ind=ind+1
+ else
+ accel(j,i)=xsolv(ind)
+ accel(j,i+nres)=xsolv(ind+1)
+ ind=ind+2
+ end if
+ enddo
+ enddo
+ enddo
+C Convert d_a to virtual-bon-vector basis
+#ifdef DEBUG
+ write (iout,*) "accel in CA-SC basis"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(accel(j,i),j=1,3),
+ & (accel(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "nnt",nnt
+#endif
+ if (nnt.eq.1) then
+ accel(:,0)=accel(:,1)
+ endif
+ do i=1,nres
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+ do j=1,3
+ accel(j,i)=accel(j,i+1)-accel(j,i)
+ enddo
+ else
+ do j=1,3
+ accel(j,i+nres)=accel(j,i+nres)-accel(j,i)
+ accel(j,i)=accel(j,i+1)-accel(j,i)
+ enddo
+ end if
+ enddo
+ accel(:,nres)=0.0d0
+ accel(:,nct)=0.0d0
+ accel(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ accel(:,0)=accel(:,1)
+ accel(:,1)=0.0d0
+ endif
+ do ichain=2,nchain
+ accel(:,chain_border1(1,ichain)-1)=
+ & accel(:,chain_border1(1,ichain))
+ accel(:,chain_border1(1,ichain))=0.0d0
+ enddo
+ do ichain=2,nchain
+ accel(:,chain_border1(1,ichain)-1)=
+ & accel(:,chain_border1(1,ichain)-1)
+ & +accel(:,chain_border(2,ichain-1))
+ accel(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+#ifdef DEBUG
+ write (iout,*) "accel in fivediaginv_mult: 1"
+ do i=0,2*nres
+ write(iout,'(i5,3f10.5)') i,(accel(j,i),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ d_a_vec(j)=accel(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_a_vec(ind+j)=accel(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ d_a_vec(ind+j)=accel(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+#ifdef TIMING
+ time_ginvmult=time_ginvmult+MPI_Wtime()-time01
+#endif
+#ifdef DEBUG
+ write (iout,*) "d_a_vec"
+ write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
+#endif
+ return
+ end
+#else