257 INTEGER(mpi),
INTENT(IN) :: n
258 REAL(mpd),
INTENT(IN OUT) :: w((n*n+n)/2)
259 REAL(mpd),
INTENT(OUT) :: aux(n)
263 aux(i)=16.0_mpd*w((i*i+i)/2)
268 IF(w(ii)+aux(i) /= aux(i))
THEN
278 w(kk+j)=w(kk+j)-w(kk+i)*ratio
305 INTEGER(mpi),
INTENT(IN) :: n
306 REAL(mpd),
INTENT(IN) :: w((n*n+n)/2)
307 REAL(mpd),
INTENT(IN) :: b(n)
308 REAL(mpd),
INTENT(OUT) :: x(n)
310 WRITE(*,*)
'before copy'
314 WRITE(*,*)
'after copy'
319 suma=suma-w(k+ii)*x(k)
324 WRITE(*,*)
'after forward'
329 suma=suma-w(kk+i)*x(k)
335 WRITE(*,*)
'after backward'
357 INTEGER(mpi),
INTENT(IN) :: n
358 REAL(mpd),
INTENT(IN) :: w((n*n+n)/2)
359 REAL(mpd),
INTENT(OUT) :: v((n*n+n)/2)
368 suma=suma-w(j+(k*k-k)/2)*v(l+(m*m-m)/2)
398 INTEGER(mpi),
INTENT(IN) :: n
399 REAL(
mpd),
INTENT(IN) :: w((n*n+n)/2)
400 REAL(
mpd),
INTENT(IN OUT) :: aux(n)
402 REAL(
mpd) :: suma,sumu,sums
407 IF(w((i*i+i)/2) < w((is*is+is)/2)) is=i
408 IF(w((i*i+i)/2) > w((ir*ir+ir)/2)) ir=i
415 IF(i == ir) suma=1.0_mpd
418 suma=suma-w(kk+i)*aux(k)
422 sumu=sumu+aux(i)*aux(i)
424 xlan=real(w((ir*ir+ir)/2)*sqrt(sumu),
mps)
425 IF(xlan /= 0.0) xlan=1.0/xlan
431 sums=sums+w(is+(i*i-i)/2)**2
435 IF(w((is*is+is)/2) /= 0.0) xla1=real(sqrt(sums)/w((is*is+is)/2),
mps)
438 IF(xla1 > 0.0.AND.xlan > 0.0) cond=xla1/xlan
469 INTEGER(mpi),
INTENT(IN) :: n
470 INTEGER(mpi),
INTENT(IN) :: mp1
471 REAL(mpd),
INTENT(IN OUT) :: w(mp1,n)
472 REAL(mpd),
INTENT(OUT) :: aux(n)
476 aux(i)=16.0_mpd*w(1,i)
479 IF(w(1,i)+aux(i) /= aux(i))
THEN
484 DO j=1,min(mp1-1,n-i)
486 DO k=1,min(mp1-1,n-i)+1-j
487 w(k,i+j)=w(k,i+j)-w(k+j,i)*rxw
511 INTEGER(mpi),
INTENT(IN) :: n
512 INTEGER(mpi),
INTENT(IN) :: mp1
513 REAL(mpd),
INTENT(IN) :: w(mp1,n)
514 REAL(mpd),
INTENT(IN) :: b(n)
515 REAL(mpd),
INTENT(OUT) :: x(n)
521 DO j=1,min(mp1-1,n-i)
522 x(j+i)=x(j+i)-w(j+1,i)*x(i)
527 DO j=1,min(mp1-1,n-i)
528 rxw=rxw-w(j+1,i)*x(j+i)
551 INTEGER(mpi),
INTENT(IN) :: n
552 INTEGER(mpi),
INTENT(IN) :: mp1
553 REAL(mpd),
INTENT(IN) :: w(mp1,n)
554 REAL(mpd),
INTENT(OUT) :: v(mp1,n)
558 DO j=i,max(1,i-mp1+1),-1
559 DO k=j+1,min(n,j+mp1-1)
560 rxw=rxw-v(1+abs(i-k),min(i,k))*w(1+k-j,j)
585 INTEGER(mpi),
INTENT(IN) :: n
586 INTEGER(mpi),
INTENT(IN) :: mp1
587 REAL(mpd),
INTENT(IN) :: w(mp1,n)
588 REAL(mpd),
INTENT(OUT) :: vs((n*n+n)/2)
592 DO j=i,max(1,i-mp1+1),-1
593 DO k=j+1,min(n,j+mp1-1)
594 rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
622 INTEGER(mpi),
INTENT(IN) :: n
623 INTEGER(mpi),
INTENT(IN) :: mp1
624 REAL(mpd),
INTENT(IN) :: w(mp1,n)
625 REAL(mpd),
INTENT(OUT) :: vs((n*n+n)/2)
630 DO k=j+1,min(n,j+mp1-1)
631 rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
657 INTEGER(mpi),
INTENT(IN) :: n
658 INTEGER(mpi),
INTENT(IN) :: mp1
659 REAL(mpd),
INTENT(IN OUT) :: w(mp1,n)
660 REAL(mpd),
INTENT(OUT) :: b(n)
663 INTEGER(mpi) :: irho(5)
672 IF(i+1-j > 0.AND.nj < 5)
THEN
674 rho=real(w(j,i+1-j)/(err*sqrt(w(1,i+1-j))),mps)
675 irho(nj)=nint(100.0*abs(rho),
mpi)
676 IF(rho < 0.0) irho(nj)=-irho(nj)
679 WRITE(*,102) i,b(i),err,(irho(j),j=1,nj)
682101
FORMAT(5x,
'i Param',7x,
'error',7x,
' c(i,i-1) c(i,i-2)'/)
683102
FORMAT(1x,i5,2g12.4,1x,5i9)
684103
FORMAT(33x,
'(correlation coefficients in percent)')
701 INTEGER(mpi),
INTENT(IN) :: mp1
702 INTEGER(mpi),
INTENT(IN) :: n
703 REAL(mpd),
INTENT(IN OUT) :: w(mp1,n)
711 WRITE(*,102) i,(w(j,i),j=1,mp1)
714101
FORMAT(5x,
'i Diag ')
715102
FORMAT(1x,i5,6g12.4)
750 INTEGER(mpi),
INTENT(IN OUT) :: n
751 REAL(mpd),
INTENT(IN OUT) :: w(2,n)
752 REAL(mpd),
INTENT(OUT) :: aux(n)
755 aux(i)=16.0_mpd*w(1,i)
758 IF(w(1,i)+aux(i) /= aux(i))
THEN
759 w(1,i)=1.0_mpd/w(1,i)
761 w(1,i+1)=w(1,i+1)-w(2,i)*rxw
768 IF(w(1,n)+aux(n) > aux(n))
THEN
769 w(1,n)=1.0_mpd/w(1,n)
790 INTEGER(mpi),
INTENT(IN) :: n
791 REAL(mpd),
INTENT(IN) :: w(2,n)
792 REAL(mpd),
INTENT(IN) :: b(n)
793 REAL(mpd),
INTENT(OUT) :: x(n)
800 x(i+1)=x(i+1)-w(2,i)*x(i)
804 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)
821 INTEGER(mpi),
INTENT(IN) :: n
822 REAL(mpd),
INTENT(IN) :: w(2,n)
823 REAL(mpd),
INTENT(OUT) :: v(2,n)
829 v(2,n-1)=-v(1,n )*w(2,n-1)
831 v(1,i )= w(1,i )-v(2,i )*w(2,i )
832 v(2,i-1)=-v(1,i )*w(2,i-1)
834 v(1,2)= w(1,2)-v(2,2)*w(2,2)
835 v(2,1)=-v(1,2)*w(2,1)
836 v(1,1)= w(1,1)-v(2,1)*w(2,1)
872 INTEGER(mpi),
INTENT(IN OUT) :: n
873 REAL(mpd),
INTENT(IN OUT) :: w(3,n)
874 REAL(mpd),
INTENT(OUT) :: aux(n)
878 aux(i)=16.0_mpd*w(1,i)
881 IF(w(1,i)+aux(i) /= aux(i))
THEN
882 w(1,i)=1.0_mpd/w(1,i)
884 w(1,i+1)=w(1,i+1)-w(2,i)*rxw
885 w(2,i+1)=w(2,i+1)-w(3,i)*rxw
888 w(1,i+2)=w(1,i+2)-w(3,i)*rxw
896 IF(w(1,n-1)+aux(n-1) > aux(n-1))
THEN
897 w(1,n-1)=1.0_mpd/w(1,n-1)
898 rxw=w(2,n-1)*w(1,n-1)
899 w(1,n)=w(1,n)-w(2,n-1)*rxw
905 IF(w(1,n)+aux(n) > aux(n))
THEN
906 w(1,n)=1.0_mpd/w(1,n)
926 INTEGER(mpi),
INTENT(IN) :: n
927 REAL(mpd),
INTENT(IN) :: w(3,n)
928 REAL(mpd),
INTENT(IN) :: b(n)
929 REAL(mpd),
INTENT(OUT) :: x(n)
935 x(i+1)=x(i+1)-w(2,i)*x(i)
936 x(i+2)=x(i+2)-w(3,i)*x(i)
938 x(n)=x(n)-w(2,n-1)*x(n-1)
940 x(n-1)=x(n-1)*w(1,n-1)-w(2,n-1)*x(n)
942 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2)
959 INTEGER(mpi),
INTENT(IN) :: n
960 REAL(mpd),
INTENT(IN) :: w(3,n)
961 REAL(mpd),
INTENT(OUT) :: v(3,n)
967 v(2,n-1)=-v(1,n )*w(2,n-1)
968 v(3,n-2)=-v(2,n-1)*w(2,n-2)-v(1,n )*w(3,n-2)
969 v(1,n-1)= w(1,n-1) -v(2,n-1)*w(2,n-1)
970 v(2,n-2)=-v(1,n-1)*w(2,n-2)-v(2,n-1)*w(3,n-2)
971 v(3,n-3)=-v(2,n-2)*w(2,n-3)-v(1,n-1)*w(3,n-3)
973 v(1,i )= w(1,i ) -v(2,i )*w(2,i )-v(3,i)*w(3,i )
974 v(2,i-1)=-v(1,i )*w(2,i-1)-v(2,i)*w(3,i-1)
975 v(3,i-2)=-v(2,i-1)*w(2,i-2)-v(1,i)*w(3,i-2)
977 v(1,2)= w(1,2) -v(2,2)*w(2,2)-v(3,2)*w(3,2)
978 v(2,1)=-v(1,2)*w(2,1)-v(2,2)*w(3,1)
979 v(1,1)= w(1,1) -v(2,1)*w(2,1)-v(3,1)*w(3,1)
1044 INTEGER(mpi),
INTENT(IN) :: n
1045 REAL(mpd),
INTENT(OUT) :: w((n*n+n)/2)
1046 INTEGER(mpi) :: i,j,k
1047 REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
1053 gamm=max(gamm,abs(w((k*k+k)/2)))
1055 xchi=max(xchi,abs(w((j*j-j)/2+k)))
1058 beta=sqrt(max(gamm,xchi/max(1.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
1059 delta=epsm*max(1.0_mpd,gamm+xchi)
1063 w((k*k-k)/2+i)=w((k*k-k)/2+i)*w((i*i+i)/2)
1067 w((j*j-j)/2+k)=w((j*j-j)/2+k)-w((k*k-k)/2+i)*w((j*j-j)/2+i)
1072 theta=max(theta,abs(w((j*j-j)/2+k)))
1074 w((k*k+k)/2)=1.0_mpd/max(abs(w((k*k+k)/2)),(theta/beta)**2,delta)
1076 w((j*j+j)/2)=w((j*j+j)/2)-w((j*j-j)/2+k)**2*w((k*k+k)/2)
1096 INTEGER(mpi),
INTENT(IN OUT) :: mp1
1097 INTEGER(mpi),
INTENT(IN) :: n
1098 REAL(mpd),
INTENT(OUT) :: w(mp1,n)
1099 INTEGER(mpi) :: i,j,k
1100 REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
1106 gamm=max(gamm,abs(w(1,k)))
1107 DO j=2,min(mp1,n-k+1)
1108 xchi=max(xchi,abs(w(j,k)))
1111 beta=sqrt(max(gamm,xchi/max(1.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
1112 delta=epsm*max(1.0_mpd,gamm+xchi)
1116 w(i,k-i+1)=w(i,k-i+1)*w(1,k-i+1)
1118 DO j=2,min(mp1,n-k+1)
1119 DO i=max(2,j+k+1-mp1),k
1120 w(j,k)=w(j,k)-w(k-i+2,i-1)*w(j-i+k+1,i-1)
1124 DO j=2,min(mp1,n-k+1)
1125 theta=max(theta,abs(w(j,k)))
1127 w(1,k)=1.0_mpd/max(abs(w(1,k)),(theta/beta)**2,delta)
1128 DO j=2,min(mp1,n-k+1)
1129 w(1,k+j-1)=w(1,k+j-1)-w(1,k)*w(j,k)**2
subroutine dbcslv(w, mp1, n, b, x)
solution B -> X
subroutine dchslv(w, n, b, x)
solution B -> X
subroutine dbcprv(w, mp1, n, b)
Print corr band matrix and pars.
subroutine dbcinv(w, mp1, n, vs)
VS = inverse symmetric matrix.
subroutine dbfdec(w, mp1, n)
Decomposition of symmetric band matrix.
subroutine dbcinb(w, mp1, n, vs)
VS = band part of inverse symmetric matrix.
real(mps) function condes(w, n, aux)
Etimate condition.
subroutine db3iel(w, n, v)
V = inverse band matrix elements.
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
subroutine dcfdec(w, n)
Decomposition of symmetric matrix.
subroutine db2dec(w, n, aux)
Decomposition (M=1).
subroutine db3dec(w, n, aux)
Decomposition (M=2).
subroutine dbcprb(w, mp1, n)
Print band matrix.
subroutine dbciel(w, mp1, n, v)
V = inverse band matrix elements.
subroutine db2slv(w, n, b, x)
solution B -> X
subroutine db3slv(w, n, b, x)
solution B -> X
subroutine dchdec(w, n, aux)
Decomposition of symmetric matrix.
subroutine dchinv(w, n, v)
inversion
subroutine db2iel(w, n, v)
V = inverse band matrix elements.
integer, parameter mpd
double precision
integer, parameter mps
single precision
integer, parameter mpi
integer