97SUBROUTINE sqminv(v,b,n,nrank,diag,next)
112 INTEGER(mpi) :: next0
114 INTEGER(mpi),
INTENT(IN) :: n
115 REAL(mpd),
INTENT(IN OUT) :: v((n*n+n)/2)
116 REAL(mpd),
INTENT(OUT) :: b(n)
117 INTEGER(mpi),
INTENT(OUT) :: nrank
118 REAL(mpd),
INTENT(OUT) :: diag(n)
119 INTEGER(mpi),
INTENT(OUT) :: next(n)
126 eps = 16.0_mpd * epsilon(eps)
132 diag(i)=abs(v((i*i+i)/2))
145 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j)))
THEN
193 v(jl)=v(jl)-v(lk)*vjk
200 IF(next(k) /= 0)
THEN
203 IF(next(j) /= 0) v((k*k-k)/2+j)=0.0_mpd
230SUBROUTINE sqminl(v,b,n,nrank,diag,next,vk,mon)
239 INTEGER(mpi) :: next0
241 INTEGER(mpi),
INTENT(IN) :: n
242 REAL(mpd),
INTENT(IN OUT) :: v((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
243 REAL(mpd),
INTENT(OUT) :: b(n)
244 INTEGER(mpi),
INTENT(OUT) :: nrank
245 REAL(mpd),
INTENT(OUT) :: diag(n)
246 INTEGER(mpi),
INTENT(OUT) :: next(n)
247 REAL(mpd),
INTENT(OUT) :: vk(n)
248 INTEGER(mpi),
INTENT(IN) :: mon
262 REAL(mpd),
PARAMETER :: eps=1.0e-10_mpd
269 diag(i)=abs(v((i8*i8+i8)/2))
284 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j)))
THEN
335 v(jl+1:jl+j)=v(jl+1:jl+j)-vk(1:j)*vjk
342 IF(next(k) /= 0)
THEN
345 IF(next(j) /= 0) v(kk+int(j,mpl))=0.0_mpd
352 10
DO jj=1,(int(n,mpl)*int(n+1,mpl))/2
374 INTEGER(mpi),
INTENT(IN) :: n
375 REAL(mpd),
INTENT(OUT) :: diag(n)
376 REAL(mpd),
INTENT(OUT) :: u(n,n)
377 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
378 REAL(mpd),
INTENT(OUT) :: work(n)
379 INTEGER(mpi),
INTENT(OUT) :: iwork(n)
382 INTEGER(mpi),
PARAMETER :: itmax=30
383 REAL(mpd),
PARAMETER :: tol=epsilon(tol)
384 REAL(mpd),
PARAMETER :: eps=epsilon(eps)
421 IF(abs(u(i,k)) > tol) g=g+u(i,k)*u(i,k)
431 IF(f >= 0.0_mpd) sh=-sh
442 IF(abs(u(j,k)) > tol.AND.abs(u(i,k)) > tol)
THEN
447 IF(abs(u(k,j)) > tol.AND.abs(u(i,k)) > tol)
THEN
462 u(j,k)=u(j,k)-f*work(k)-g*u(i,k)
474 IF(diag(i) /= 0.0)
THEN
481 u(k,j)=u(k,j)-g*u(k,i)
503 h=eps*(abs(diag(l))+abs(work(l)))
506 IF(abs(work(m)) <= b)
GO TO 10
50910
IF(m == l)
GO TO 30
51120
IF(j == itmax)
THEN
512 WRITE(*,*)
'DEVROT: Iteration limit reached'
513 CALL peend(32,
'Aborted, iteration limit reached in diagonalization')
518 p=(diag(l+1)-g)/(2.0_mpd*work(l))
521 IF(p < 0.0_mpd) diag(l)=diag(l)/(p-r)
522 IF(p >= 0.0_mpd) diag(l)=diag(l)/(p+r)
535 IF(abs(p) >= abs(work(i)))
THEN
544 work(i+1)=s*work(i)*r
549 diag(i+1)=h+s*(c*g+s*diag(i))
553 u(k,i+1)=s*u(k,i)+c*h
559 IF(abs(work(l)) > b)
GO TO 20
57260
IF(diag(iwork(l+m)) > diag(iwork(l)))
THEN
583 IF(iwork(i) /= i)
THEN
616 INTEGER(mpi),
INTENT(IN) :: n
617 REAL(mpd),
INTENT(IN) :: diag(n)
618 REAL(mpd),
INTENT(IN) :: u(n,n)
619 REAL(mpd),
INTENT(IN) :: b(n)
620 REAL(mpd),
INTENT(OUT) :: coef(n)
627 IF(diag(i) > 0.0_mpd)
THEN
630 dsum=dsum+u(j,i)*b(j)
632 coef(i)=abs(dsum)/sqrt(diag(i))
654 INTEGER(mpi),
INTENT(IN) :: n
655 REAL(mpd),
INTENT(IN) :: diag(n)
656 REAL(mpd),
INTENT(IN) :: u(n,n)
657 REAL(mpd),
INTENT(IN) :: b(n)
658 REAL(mpd),
INTENT(OUT) :: x(n)
659 REAL(mpd),
INTENT(OUT) :: work(n)
668 IF(diag(j) /= 0.0_mpd)
THEN
705 INTEGER(mpi),
INTENT(IN) :: n
706 REAL(mpd),
INTENT(IN) :: diag(n)
707 REAL(mpd),
INTENT(IN) :: u(n,n)
708 REAL(mpd),
INTENT(OUT) :: v((n*n+n)/2)
717 IF(diag(k) /= 0.0_mpd)
THEN
718 dsum=dsum+u(i,k)*u(j,k)/diag(k)
756 INTEGER(mpi),
INTENT(IN) :: n
757 REAL(mpd),
INTENT(IN OUT) :: g((n*n+n)/2)
764 IF(g(ii) /= 0.0) g(ii)=1.0/g(ii)
770 g(kk+j)=g(kk+j)-g(kk+i)*ratio
801 INTEGER(mpi),
INTENT(IN) :: n
802 REAL(mpd),
INTENT(IN) :: g((n*n+n)/2)
803 REAL(mpd),
INTENT(IN OUT) :: x(n)
809 dsum=dsum-g(k+ii)*x(k)
818 dsum=dsum-g(kk+i)*x(k)
849 INTEGER(mpi),
INTENT(IN) :: n
850 REAL(mpd),
INTENT(IN) :: g((n*n+n)/2)
851 REAL(mpd),
INTENT( OUT) :: v((n*n+n)/2)
860 dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2)
891SUBROUTINE chdec2(g,n,nrank,evmax,evmin,mon)
902 INTEGER(mpi),
INTENT(IN) :: n
903 REAL(mpd),
INTENT(IN OUT) :: g((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
904 INTEGER(mpi),
INTENT(OUT) :: nrank
905 REAL(mpd),
INTENT(OUT) :: evmin
906 REAL(mpd),
INTENT(OUT) :: evmax
907 INTEGER(mpi),
INTENT(IN) :: mon
910 ii=(int(n,mpl)*int(n+1,mpl))/2
913 IF(mon>0)
CALL monpgs(n+1-i)
914 IF (g(ii) > 0.0_mpd)
THEN
921 evmax=max(evmax,g(ii))
922 evmin=min(evmin,g(ii))
933 ratio=g(ii+j)*g(ii+i)
934 IF (ratio == 0.0_mpd) cycle
935 jj=(int(j-1,mpl)*int(j,mpl))/2
936 g(jj+1:jj+j)=g(jj+1:jj+j)-g(ii+1:ii+j)*ratio
939 g(ii+1:ii+i-1)=g(ii+1:ii+i-1)*g(ii+i)
963 INTEGER(mpi),
INTENT(IN) :: n
964 REAL(mpd),
INTENT(IN) :: g((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
965 REAL(mpd),
INTENT(IN OUT) :: x(n)
967 ii=(int(n,mpl)*int(n+1,mpl))/2
972 dsum=dsum-g(kk+i)*x(k)
981 dsum=dsum-g(k+ii)*x(k)
1029 INTEGER(mpi),
INTENT(IN) :: n
1030 REAL(mpd),
INTENT(IN OUT) :: val(*)
1031 INTEGER(mpi),
INTENT(IN) :: ilptr(n)
1037 REAL(mpd),
PARAMETER :: one=1.0_mpd
1038 REAL(mpd),
PARAMETER :: two=2.0_mpd
1039 REAL(mpd),
PARAMETER :: eps = epsilon(eps)
1041 WRITE(*,*)
'Variable band matrix Cholesky decomposition'
1046 IF(ilptr(i) == j)
THEN
1047 IF(val(j) <= 0.0_mpd)
GO TO 01
1048 dgamma=max(dgamma,abs(val(j)))
1054 WRITE(*,*)
' ',in,
' positive diagonal elements'
1059 IF(ilptr(i) == j)
THEN
1062 xi=max(xi,abs(val(j)))
1066 delta=eps*max(1.0_mpd,dgamma+xi)
1068 IF(n > 1) sn=1.0_mpd/sqrt(real(n*n-1,mpd))
1069 beta=sqrt(max(eps,dgamma,xi*sn))
1070 WRITE(*,*)
' DELTA and BETA ',delta,beta
1073 mk=k-ilptr(k)+ilptr(k-1)+1
1078 mj=j-ilptr(j)+ilptr(j-1)+1
1083 -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i)
1087 theta=max(theta,abs(val(kj)))
1090 IF(val(ilptr(j)) /= 0.0_mpd)
THEN
1091 val(kj)=val(kj)/val(ilptr(j))
1100 val(kj)=max(abs(val(kj)),(theta/beta)**2,delta)
1101 IF(valkj /= val(kj))
THEN
1102 WRITE(*,*)
' Index K=',k
1103 WRITE(*,*)
' ',valkj,val(kj), (theta/beta)**2,delta,theta
1112 IF(val(ilptr(k)) /= 0.0_mpd) val(ilptr(k))=1.0_mpd/val(ilptr(k))
1132 INTEGER(mpi),
INTENT(IN) :: n
1133 REAL(mpd),
INTENT(IN OUT) :: val(*)
1134 INTEGER(mpi),
INTENT(IN) :: ilptr(n)
1139 IF(val(ilptr(k)) > val(ilptr(ks))) ks=k
1140 IF(val(ilptr(k)) < val(ilptr(kr))) kr=k
1141 IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k
1143 WRITE(*,*)
' Index value ',ks,val(ilptr(ks))
1144 WRITE(*,*)
' Index value ',kp,val(ilptr(kp))
1145 WRITE(*,*)
' Index value ',kr,val(ilptr(kr))
1169 INTEGER(mpi),
INTENT(IN) :: n
1170 REAL(mpd),
INTENT(IN OUT) :: val(*)
1171 INTEGER(mpi),
INTENT(IN) :: ilptr(n)
1172 REAL(mpd),
INTENT(IN OUT) :: x(n)
1175 mk=k-ilptr(k)+ilptr(k-1)+1
1177 x(k)=x(k)-val(ilptr(k)-k+j)*x(j)
1182 x(k)=x(k)*val(ilptr(k))
1186 mk=k-ilptr(k)+ilptr(k-1)+1
1188 x(j)=x(j)-val(ilptr(k)-k+j)*x(k)
1209 INTEGER(mpi),
INTENT(IN) :: n
1210 REAL(
mpd),
INTENT(IN) :: dx(n)
1211 REAL(
mpd),
INTENT(IN) :: dy(n)
1215 dtemp=dtemp+dx(i)*dy(i)
1217 DO i =mod(n,5)+1,n,5
1218 dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) &
1219 +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4)
1239 INTEGER(mpi),
INTENT(IN) :: n
1240 REAL(mpd),
INTENT(IN) :: dx(n)
1241 REAL(mpd),
INTENT(IN OUT) :: dy(n)
1242 REAL(mpd),
INTENT(IN) :: da
1245 dy(i)=dy(i)+da*dx(i)
1248 dy(i )=dy(i )+da*dx(i )
1249 dy(i+1)=dy(i+1)+da*dx(i+1)
1250 dy(i+2)=dy(i+2)+da*dx(i+2)
1251 dy(i+3)=dy(i+3)+da*dx(i+3)
1276 INTEGER(mpi),
INTENT(IN) :: n
1277 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
1278 REAL(mpd),
INTENT(IN) :: a(n)
1279 REAL(mpd),
INTENT(OUT) :: b(n)
1287 dsum=dsum+v(ij)*a(j)
1318 INTEGER(mpi),
INTENT(IN) :: n
1319 REAL(mpd),
INTENT(IN) :: v((INT(n,mpl)*INT(n,mpl)+INT(n,mpl))/2)
1320 REAL(mpd),
INTENT(IN) :: a(n)
1321 REAL(mpd),
INTENT(OUT) :: b(n)
1331 dsum=dsum+v(ij)*a(j)
1359 INTEGER(mpi),
INTENT(IN) :: m
1360 INTEGER(mpi),
INTENT(IN) :: n
1361 REAL(mpd),
INTENT(IN) :: a(n*m)
1362 REAL(mpd),
INTENT(IN) :: x(n)
1363 REAL(mpd),
INTENT(OUT) :: y(m)
1371 y(i)=y(i)+a(ij)*x(j)
1404 INTEGER(mpi),
INTENT(IN) :: n
1405 INTEGER(mpi),
INTENT(IN) :: m
1406 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
1407 REAL(mpd),
INTENT(IN) :: a(n*m)
1408 REAL(mpd),
INTENT(INOUT) :: w((m*m+m)/2)
1410 INTEGER(mpi),
INTENT(IN) :: iopt
1432 cik=cik+a(il+l)*v(lk)
1436 cik=cik+a(il+l)*v(lk)
1442 w(ij)=w(ij)+cik*a(jk)
1485 INTEGER(mpi),
INTENT(IN) :: n
1486 INTEGER(mpi),
INTENT(IN) :: m
1487 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
1488 REAL(mpd),
INTENT(IN) :: a(n*m)
1489 REAL(mpd),
INTENT(INOUT) :: w((m*m+m)/2)
1490 INTEGER(mpi),
INTENT(IN) :: is(2*n*m+n+m+1)
1491 INTEGER(mpi),
INTENT(IN) :: iopt
1492 INTEGER(mpi),
INTENT(OUT) :: sc(n)
1513 DO l=is(i)+1,is(i+1),2
1516 lk=sc(max(k,ic))+min(k,ic)
1519 DO j=is(m+k)+1,is(m+k+1),2
1524 w(ij)=w(ij)+cik*a(in)
1551 INTEGER(mpi) :: mc(15)
1556 INTEGER(mpi),
INTENT(IN) :: lun
1557 INTEGER(mpi),
INTENT(IN) :: n
1558 REAL(mpd),
INTENT(IN) :: x(n)
1559 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
1568 IF(v(ii) > 0.0) err=sqrt(real(v(ii),mps))
1575 pd=real(v(ii)*v(jj),mps)
1576 IF(pd > 0.0) rho=real(v(ij),mps)/sqrt(pd)
1578 mc(l)=nint(100.0*abs(rho),
mpi)
1579 IF(rho < 0.0) mc(l)=-mc(l)
1580 IF(j == i.OR.l == 15)
THEN
1583 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l-1)
1585 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l)
1589 WRITE(lun,103) (mc(m),m=1,l-1)
1591 WRITE(lun,103) (mc(m),m=1,l)
1601101
FORMAT(9x,
'Param',7x,
'error',7x,
'correlation coefficients'/)
1602102
FORMAT(1x,i5,2g12.4,1x,15i5)
1604104
FORMAT(33x,
'(correlation coefficients in percent)')
1619 INTEGER(mpi),
PARAMETER :: istp=6
1627 INTEGER(mpi),
INTENT(IN) :: lun
1628 INTEGER(mpi),
INTENT(IN) :: n
1629 REAL(mpd),
INTENT(IN) :: v((n*n+n)/2)
1639 WRITE(lun,102) i, ip+1-ips, (v(k),k=ip+1,min(ipn,ipe))
1646101
FORMAT(1x,
'--- DBPRV -----------------------------------')
1647102
FORMAT(1x,2i3,6g12.4)
1670 INTEGER(mpi),
INTENT(IN) :: n
1671 REAL(mps),
INTENT(IN OUT) :: a(n)
1692 IF(a(j) < a(j+1)) j=j+1
1718 INTEGER(mpi) :: nlev
1719 parameter(nlev=2*32)
1725 INTEGER(mpi) :: lr(nlev)
1727 INTEGER(mpi) :: maxlev
1731 INTEGER(mpi),
INTENT(IN) :: n
1732 INTEGER(mpi),
INTENT(IN OUT) :: a(n)
1740 IF(a(l) > a(r))
THEN
1759 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1764 IF(a(i) < a1)
GO TO 20
1766 IF(a(j) > a1)
GO TO 30
1773 IF(lev+2 > nlev)
THEN
1774 CALL peend(33,
'Aborted, stack overflow in quicksort')
1775 stop
'SORT1K (quicksort): stack overflow'
1787 maxlev=max(maxlev,lev)
1803 INTEGER(mpi) :: nlev
1804 parameter(nlev=2*32)
1810 INTEGER(mpi) ::lr(nlev)
1812 INTEGER(mpi) ::maxlev
1817 INTEGER(mpi),
INTENT(IN OUT) :: a(2,*)
1818 INTEGER(mpi),
INTENT(IN) :: n
1825 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r)))
THEN
1837 WRITE(*,*)
'SORT2K (quicksort): maxlevel used/available =', maxlev,
'/64'
1846 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1852 IF(a(1,i) < a1)
GO TO 20
1853 IF(a(1,i) == a1.AND.a(2,i) < a2)
GO TO 20
1855 IF(a(1,j) > a1)
GO TO 30
1856 IF(a(1,j) == a1.AND.a(2,j) > a2)
GO TO 30
1866 IF(lev+2 > nlev)
THEN
1867 CALL peend(33,
'Aborted, stack overflow in quicksort')
1868 stop
'SORT2K (quicksort): stack overflow'
1880 maxlev=max(maxlev,lev)
1896 INTEGER(mpi) :: nlev
1897 parameter(nlev=2*32)
1903 INTEGER(mpi) ::lr(nlev)
1905 INTEGER(mpi) ::maxlev
1908 INTEGER(mpi) ::at(3)
1910 INTEGER(mpi),
INTENT(IN OUT) :: a(3,*)
1911 INTEGER(mpi),
INTENT(IN) :: n
1918 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r)))
THEN
1927 WRITE(*,*)
'SORT2I (quicksort): maxlevel used/available =', maxlev,
'/64'
1936 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1942 IF(a(1,i) < a1)
GO TO 20
1943 IF(a(1,i) == a1.AND.a(2,i) < a2)
GO TO 20
1945 IF(a(1,j) > a1)
GO TO 30
1946 IF(a(1,j) == a1.AND.a(2,j) > a2)
GO TO 30
1948 IF(a(1,i) == a(1,j).AND.a(2,i) == a(2,j))
GO TO 20
1954 IF(lev+2 > nlev)
THEN
1955 CALL peend(33,
'Aborted, stack overflow in quicksort')
1956 stop
'SORT2I (quicksort): stack overflow'
1968 maxlev=max(maxlev,lev)
1985 INTEGER(mpi) :: nlev
1986 parameter(nlev=2*32)
1992 INTEGER(mpi) ::lr(nlev)
1994 INTEGER(mpi) ::maxlev
1997 INTEGER(mpi) ::at(4)
2000 INTEGER(mpi),
INTENT(IN OUT) :: a(4,*)
2001 INTEGER(mpl),
INTENT(IN OUT) :: b(*)
2002 INTEGER(mpi),
INTENT(IN) :: n
2010 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r)))
THEN
2022 WRITE(*,*)
'SORT22l (quicksort): maxlevel used/available =', maxlev,
'/64'
2031 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
2037 IF(a(1,i) < a1)
GO TO 20
2038 IF(a(1,i) == a1.AND.a(2,i) < a2)
GO TO 20
2040 IF(a(1,j) > a1)
GO TO 30
2041 IF(a(1,j) == a1.AND.a(2,j) > a2)
GO TO 30
2051 IF(lev+2 > nlev)
THEN
2052 CALL peend(33,
'Aborted, stack overflow in quicksort')
2053 stop
'SORT22l (quicksort): stack overflow'
2065 maxlev=max(maxlev,lev)
2083 INTEGER(mpi),
INTENT(IN) :: n
2084 INTEGER(mpi),
INTENT(IN) :: nd
2087 REAL(
mps) ::table(30,3)
2090 DATA sn/0.47523,1.690140,2.782170/
2091 DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
2092 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
2093 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
2094 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
2095 4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
2096 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
2097 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
2098 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
2099 9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
2100 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
2101 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
2102 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
2112 chindl=(sn(m)+sqrt(real(nd+nd-1,
mps)))**2/real(nd+nd,
mps)
2167 INTEGER(mpi),
INTENT(IN) :: n
2168 INTEGER(mpi),
INTENT(IN) :: india(n)
2169 INTEGER(mpi),
INTENT(OUT) :: nrkd
2170 INTEGER(mpi),
INTENT(IN) :: iopt
2171 REAL(mpd),
INTENT(IN OUT) :: c(india(n))
2174 eps = 16.0_mpd * epsilon(eps)
2179 IF(c(india(1)) > 0.0)
THEN
2180 c(india(1))=1.0_mpd/sqrt(c(india(1)))
2187 mk=k-india(k)+india(k-1)+1
2190 IF(j > 1) mj=j-india(j)+india(j-1)+1
2196 c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i)
2199 IF(j /= k) c(kj)=c(kj)*diag
2202 IF(c(india(k)) > eps*diag)
THEN
2203 c(india(k))=1.0_mpd/sqrt(c(india(k)))
2206 c(india(k)-k+j)=0.0_mpd
2208 IF (iopt > 0 .and. diag > 0.0)
THEN
2209 c(india(k))=1.0_mpd/sqrt(diag)
2244 INTEGER(mpi),
INTENT(IN) :: n
2245 INTEGER(mpi),
INTENT(IN) :: india(n)
2246 REAL(mpd),
INTENT(IN OUT) :: c(india(n))
2247 REAL(mpd),
INTENT(IN OUT) :: x(n)
2249 x(1)=x(1)*c(india(1))
2251 DO j=k-india(k)+india(k-1)+1,k-1
2252 x(k)=x(k)-c(india(k)-k+j)*x(j)
2254 x(k)=x(k)*c(india(k))
2283 INTEGER(mpi),
INTENT(IN) :: n
2284 INTEGER(mpi),
INTENT(IN) :: india(n)
2285 REAL(mpd),
INTENT(IN OUT) :: c(india(n))
2286 REAL(mpd),
INTENT(IN OUT) :: x(n)
2287 INTEGER(mpi),
INTENT(IN) :: i0
2288 INTEGER(mpi),
INTENT(IN) :: ns
2293 DO j=max(1,k-india(l)+india(l-1)+1),k-1
2294 x(k)=x(k)-c(india(l)-k+j)*x(j)
2297 x(k)=x(k)*c(india(l))
2322 INTEGER(mpi),
INTENT(IN) :: n
2323 INTEGER(mpi),
INTENT(IN) :: india(n)
2324 REAL(mpd),
INTENT(IN OUT) :: c(india(n))
2325 REAL(mpd),
INTENT(IN OUT) :: x(n)
2328 x(k)=x(k)*c(india(k))
2329 DO j=k-india(k)+india(k-1)+1,k-1
2330 x(j)=x(j)-c(india(k)-k+j)*x(k)
2333 x(1)=x(1)*c(india(1))
2365 INTEGER(mpi),
INTENT(IN) :: n
2366 INTEGER(mpi),
INTENT(IN) :: m
2367 INTEGER(mpi),
INTENT(IN) :: ls
2368 INTEGER(mpi),
INTENT(IN OUT) :: india(n+m)
2369 INTEGER(mpi),
INTENT(OUT) :: nrkd
2370 INTEGER(mpi),
INTENT(OUT) :: nrkd2
2371 REAL(mpd),
INTENT(IN OUT) :: c(india(n)+n*m+(m*m+m)/2)
2379 CALL lltdec(n,c,india,nrkd,ls)
2383 CALL lltfwd(n,c,india,c(india(n)+int(i-1,
mpl)*nn+1))
2386 jk=india(n)+nn*int(m,
mpl)
2392 c(jk)=c(jk)+c(india(n)+int(j-1,
mpl)*nn+i)*c(india(n)+int(k-1,
mpl)*nn+i)
2399 india(n+i)=india(n+i-1)+min(i,m)
2401 CALL lltdec(m,c(india(n)+nn*int(m,
mpl)+1),india(n+1),nrkd2,0)
2430 INTEGER(mpi),
INTENT(IN) :: n
2431 INTEGER(mpi),
INTENT(IN) :: m
2432 INTEGER(mpi),
INTENT(IN) :: india(n+m)
2433 REAL(mpd),
INTENT(IN OUT) :: c(india(n)+n*m+(m*m+m)/2)
2434 REAL(mpd),
INTENT(IN OUT) :: x(n+m)
2442 x(n+i)=x(n+i)-x(j)*c(india(n)+int(i-1,
mpl)*nn+j)
2445 CALL lltfwd(m,c(india(n)+nn*int(m,
mpl)+1),india(n+1),x(n+1))
2448 CALL lltbwd(m,c(india(n)+nn*int(m,
mpl)+1),india(n+1),x(n+1))
2455 x(i)=x(i)-x(n+j)*c(india(n)+int(j-1,
mpl)*nn+i)
2486SUBROUTINE equdecs(n,m,b,nm,ls,c,india,l,nrkd,nrkd2)
2510 INTEGER(mpi),
INTENT(IN) :: n
2511 INTEGER(mpi),
INTENT(IN) :: m
2512 INTEGER(mpi),
INTENT(IN) :: b
2513 INTEGER(mpl),
INTENT(IN) :: nm
2514 INTEGER(mpi),
INTENT(IN) :: ls
2515 INTEGER(mpi),
INTENT(IN OUT) :: india(n+m)
2516 INTEGER(mpi),
INTENT(IN) :: l(4*b)
2517 INTEGER(mpi),
INTENT(OUT) :: nrkd
2518 INTEGER(mpi),
INTENT(OUT) :: nrkd2
2519 REAL(mpd),
INTENT(IN OUT) :: c(nm+india(n)+(m*m+m)/2)
2526 CALL lltdec(n,c,india,nrkd,ls)
2529 ij=india(n)+(m*m+m)/2
2536 CALL lltfwds(n,c,india,c(ij+1),i0,in)
2541 ij=india(n)+(m*m+m)/2
2550 jk=india(n)+(j*j-j)/2+p1
2552 c(jk)=c(jk)+c(ij+(j-p1)*in)*c(ij+(k-p1)*in)
2571 jk=india(n)+(j*j-j)/2+p1
2573 c(jk)=c(jk)+c(jj+(j-q1)*jn+i)*c(kk+(k-p1)*in)
2585 india(n+i)=india(n+i-1)+min(i,m)
2588 CALL lltdec(m,c(india(n)+1),india(n+1),nrkd2,0)
2627 INTEGER(mpi),
INTENT(IN) :: n
2628 INTEGER(mpi),
INTENT(IN) :: m
2629 INTEGER(mpi),
INTENT(IN) :: b
2630 INTEGER(mpl),
INTENT(IN) :: nm
2631 INTEGER(mpi),
INTENT(IN) :: india(n+m)
2632 INTEGER(mpi),
INTENT(IN) :: l(4*b)
2633 REAL(mpd),
INTENT(IN OUT) :: c(nm+india(n)+(m*m+m)/2)
2634 REAL(mpd),
INTENT(IN OUT) :: x(n+m)
2639 ij=india(n)+(m*m+m)/2
2648 x(n+j)=x(n+j)-c(ij)*x(i+i0)
2653 CALL lltfwd(m,c(india(n)+1),india(n+1),x(n+1))
2655 CALL lltbwd(m,c(india(n)+1),india(n+1),x(n+1))
2660 ij=india(n)+(m*m+m)/2
2669 x(i+i0)=x(i+i0)-c(ij)*x(n+j)
2722 INTEGER(mpi),
INTENT(IN) :: p
2723 INTEGER(mpi),
INTENT(IN) :: n
2724 REAL(mpd),
INTENT(IN) :: c(n)
2725 REAL(mpd),
INTENT(OUT) :: cu(n)
2726 REAL(mpd),
INTENT(IN OUT) :: a(n,p)
2727 REAL(mpd),
INTENT(OUT) :: s((p*p+p)/2)
2728 INTEGER(mpi),
INTENT(OUT) :: nrkd
2740 IF (div > 0.0_mpd)
THEN
2741 cu(i)=1.0_mpd/sqrt(div)
2750 s(jk)=s(jk)+a(i,j)*a(i,k)
2758 IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2764 s(kk+j)=s(kk+j)-s(kk+i)*ratio
2794 INTEGER(mpi),
INTENT(IN) :: p
2795 INTEGER(mpi),
INTENT(IN) :: n
2797 REAL(mpd),
INTENT(IN) :: cu(n)
2798 REAL(mpd),
INTENT(IN) :: a(n,p)
2799 REAL(mpd),
INTENT(IN) :: s((p*p+p)/2)
2800 REAL(mpd),
INTENT(OUT) :: x(n+p)
2801 REAL(mpd),
INTENT(IN) :: y(n+p)
2811 x(n+j)=x(n+j)-a(i,j)*x(i)
2819 dsum=dsum-s(k+jj)*x(n+k)
2829 dsum=dsum+s(kk+j)*x(n+k)
2838 x(i)=x(i)-a(i,j)*x(n+j)
2899 INTEGER(mpi),
INTENT(IN) :: p
2900 INTEGER(mpi),
INTENT(IN) :: n
2901 INTEGER(mpi),
INTENT(IN) :: b
2902 INTEGER(mpl),
INTENT(IN) :: nm
2903 REAL(mpd),
INTENT(IN) :: c(n)
2904 REAL(mpd),
INTENT(OUT) :: cu(n)
2905 REAL(mpd),
INTENT(IN OUT) :: a(nm)
2906 INTEGER(mpi),
INTENT(IN) :: l(p)
2907 REAL(mpd),
INTENT(OUT) :: s((p*p+p)/2)
2908 INTEGER(mpi),
INTENT(OUT) :: nrkd
2919 IF (div > 0.0_mpd)
THEN
2920 cu(i)=1.0_mpd/sqrt(div)
2937 a(ij+(j-p1)*np)=a(ij+(j-p1)*np)*cu(i+i0)
2939 s(jk)=s(jk)+a(ij+(j-p1)*np)*a(ij+(k-p1)*np)
2950 IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2956 s(kk+j)=s(kk+j)-s(kk+i)*ratio
2996 INTEGER(mpi),
INTENT(IN) :: p
2997 INTEGER(mpi),
INTENT(IN) :: n
2998 INTEGER(mpi),
INTENT(IN) :: b
2999 INTEGER(mpl),
INTENT(IN) :: nm
3001 REAL(mpd),
INTENT(IN) :: cu(n)
3002 REAL(mpd),
INTENT(IN) :: a(nm)
3003 INTEGER(mpi),
INTENT(IN) :: l(p)
3004 REAL(mpd),
INTENT(IN) :: s((p*p+p)/2)
3005 REAL(mpd),
INTENT(OUT) :: x(n+p)
3006 REAL(mpd),
INTENT(IN) :: y(n+p)
3025 x(n+j)=x(n+j)-a(ij)*x(i+i0)
3034 dsum=dsum-s(k+jj)*x(n+k)
3044 dsum=dsum+s(kk+j)*x(n+k)
3060 x(i+i0)=x(i+i0)-a(ij)*x(n+j)
3118SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag,evdmin,evdmax)
3132 INTEGER(mpi) :: ioff
3140 INTEGER(mpi) :: joff
3144 INTEGER(mpi) :: npri
3145 INTEGER(mpi) :: nrankb
3147 INTEGER(mpi),
INTENT(IN) :: n
3148 REAL(mpd),
INTENT(IN OUT) :: v((n*n+n)/2)
3149 REAL(mpd),
INTENT(OUT) :: b(n)
3150 INTEGER(mpi),
INTENT(IN) :: nbdr
3151 INTEGER(mpi),
INTENT(IN) :: nbnd
3152 INTEGER(mpi),
INTENT(IN) :: inv
3153 INTEGER(mpi),
INTENT(OUT) :: nrank
3155 REAL(mpd),
INTENT(OUT) :: vbnd(n*(nbnd+1))
3156 REAL(mpd),
INTENT(OUT) :: vbdr(n*nbdr)
3157 REAL(mpd),
INTENT(OUT) :: aux(n*nbdr)
3158 REAL(mpd),
INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
3159 REAL(mpd),
INTENT(OUT) :: vzru(nbdr)
3160 REAL(mpd),
INTENT(OUT) :: scdiag(nbdr)
3161 INTEGER(mpi),
INTENT(OUT) :: scflag(nbdr)
3162 REAL(mpd),
INTENT(OUT) :: evdmin
3163 REAL(mpd),
INTENT(OUT) :: evdmax
3176 DO j=i,min(n,i+nbnd)
3180 vbnd(ib+(i-nb1)*mp1)=v(ip)
3198 CALL dbcdec(vbnd,mp1,nmb,aux)
3206 IF (vbnd(ip) <= 0.0_mpd)
THEN
3209 IF (vbnd(ip) == 0.0_mpd)
THEN
3210 print *,
' SQMIBB matrix singular', n, nbdr, nbnd
3212 print *,
' SQMIBB matrix not positive definite', n, nbdr, nbnd
3226 IF (nrank == 1)
THEN
3230 evdmin=min(evdmin,vbnd(ip))
3231 evdmax=max(evdmax,vbnd(ip))
3239 CALL dbcslv(vbnd,mp1,nmb,b,b)
3242 CALL dbcinv(vbnd,mp1,nmb,v)
3244 CALL dbcinb(vbnd,mp1,nmb,v)
3253 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff))
3257 vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i)
3262 CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1))
3272 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
3279 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
3280 IF (nrankb == nbdr)
THEN
3284 IF (npri >= 0) print *,
' SQMIBB undef border ', n, nbdr, nbnd, nrankb
3288 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3297 b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i)
3304 CALL dbcinv(vbnd,mp1,nmb,v)
3306 CALL dbcinb(vbnd,mp1,nmb,v)
3313 IF (inv == 1) j0=max(0,i-nbnd)
3321 ij=(ij*ij-ij)/2+min(ib,jb)
3322 v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j)
3338 ij=(ij*ij-ij)/2+min(ib,jb)
3339 v(ip2)=v(ip2)-vbk(ij)*aux(i+joff)
3346 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3390SUBROUTINE sqmibb2(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag,evdmin,evdmax)
3404 INTEGER(mpi) :: ioff
3411 INTEGER(mpi) :: joff
3412 INTEGER(mpi) :: koff
3416 INTEGER(mpi) :: npri
3417 INTEGER(mpi) :: nrankb
3419 INTEGER(mpi),
INTENT(IN) :: n
3420 REAL(mpd),
INTENT(IN OUT) :: v((n*n+n)/2)
3421 REAL(mpd),
INTENT(OUT) :: b(n)
3422 INTEGER(mpi),
INTENT(IN) :: nbdr
3423 INTEGER(mpi),
INTENT(IN) :: nbnd
3424 INTEGER(mpi),
INTENT(IN) :: inv
3425 INTEGER(mpi),
INTENT(OUT) :: nrank
3427 REAL(mpd),
INTENT(OUT) :: vbnd(n*(nbnd+1))
3428 REAL(mpd),
INTENT(OUT) :: vbdr(n*nbdr)
3429 REAL(mpd),
INTENT(OUT) :: aux(n*nbdr)
3430 REAL(mpd),
INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
3431 REAL(mpd),
INTENT(OUT) :: vzru(nbdr)
3432 REAL(mpd),
INTENT(OUT) :: scdiag(nbdr)
3433 INTEGER(mpi),
INTENT(OUT) :: scflag(nbdr)
3434 REAL(mpd),
INTENT(OUT) :: evdmin
3435 REAL(mpd),
INTENT(OUT) :: evdmax
3448 DO j=i,min(nmb,i+nbnd)
3452 vbnd(ib+(i-1)*mp1)=v(ip)
3461 vbdr(ioff+j)=v(ip+j)
3467 CALL dbcdec(vbnd,mp1,nmb,aux)
3475 IF (vbnd(ip) <= 0.0_mpd)
THEN
3478 IF (vbnd(ip) == 0.0_mpd)
THEN
3479 print *,
' SQMIBB2 matrix singular', n, nbdr, nbnd
3481 print *,
' SQMIBB2 matrix not positive definite', n, nbdr, nbnd
3495 IF (nrank == 1)
THEN
3499 evdmin=min(evdmin,vbnd(ip))
3500 evdmax=max(evdmax,vbnd(ip))
3507 CALL dbcslv(vbnd,mp1,nmb,b,b)
3510 CALL dbcinv(vbnd,mp1,nmb,v)
3512 CALL dbcinb(vbnd,mp1,nmb,v)
3521 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff+1),aux(ioff+1))
3525 vzru(ib)=vzru(ib)-b(i)*aux(ioff+i)
3530 CALL dbcslv(vbnd,mp1,nmb,b,b)
3539 vbk(ip)=vbdr(koff+jb)
3541 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
3550 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
3551 IF (nrankb == nbdr)
THEN
3555 IF (npri >= 0) print *,
' SQMIBB2 undef border ', n, nbdr, nbnd, nrankb
3559 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
3567 b(i)=b(i)-vzru(ib)*aux(ioff+i)
3575 CALL dbcinv(vbnd,mp1,nmb,v)
3577 CALL dbcinb(vbnd,mp1,nmb,v)
3585 IF (inv == 1) j0=max(0,i-nbnd)
3592 ij=(ij*ij-ij)/2+min(ib,jb)
3593 v(ip1)=v(ip1)+vbk(ij)*aux(ioff+i)*aux(joff+j)
3612 ij=(ij*ij-ij)/2+min(ib,jb)
3613 v(ip1)=v(ip1)-vbk(ij)*aux(i+joff)
subroutine dbcslv(w, mp1, n, b, x)
solution B -> X
subroutine dbcinv(w, mp1, n, vs)
VS = inverse symmetric matrix.
subroutine dbcinb(w, mp1, n, vs)
VS = band part of inverse symmetric matrix.
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
subroutine monpgs(i)
Progress monitoring.
subroutine cholin(g, v, n)
Inversion after decomposition.
subroutine vabmmm(n, val, ilptr)
Variable band matrix print minimum and maximum.
subroutine dbavat(v, a, w, n, m, iopt)
A V AT product (similarity).
subroutine sqminl(v, b, n, nrank, diag, next, vk, mon)
Matrix inversion for LARGE matrices.
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
subroutine sort22l(a, b, n)
Quick sort 2 with index.
subroutine dbmprv(lun, x, v, n)
Print symmetric matrix, vector.
subroutine dbaxpy(n, da, dx, dy)
Multiply, addition.
subroutine dbavats(v, a, is, w, n, m, iopt, sc)
A V AT product (similarity, sparse).
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
Bordered band matrix.
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
subroutine heapf(a, n)
Heap sort direct (real).
real(mpd) function dbdot(n, dx, dy)
Dot product.
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
subroutine vabdec(n, val, ilptr)
Variable band matrix decomposition.
subroutine chslv2(g, x, n)
Solve A*x=b using Cholesky decomposition.
subroutine vabslv(n, val, ilptr, x)
Variable band matrix solution.
subroutine choldc(g, n)
Cholesky decomposition.
subroutine sort1k(a, n)
Quick sort 1.
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
subroutine lltbwd(n, c, india, x)
Backward solution.
subroutine presols(p, n, b, nm, cu, a, l, s, x, y)
Constrained (sparse) preconditioner, solution.
subroutine dbgax(a, x, y, m, n)
Multiply general M-by-N matrix A and N-vector X.
subroutine lltfwds(n, c, india, x, i0, ns)
Forward solution (sparse).
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
subroutine dbprv(lun, v, n)
Print symmetric matrix.
subroutine devinv(n, diag, u, v)
Inversion by diagonalization.
subroutine lltdec(n, c, india, nrkd, iopt)
LLT decomposition.
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
Band bordered matrix.
subroutine equdecs(n, m, b, nm, ls, c, india, l, nrkd, nrkd2)
Decomposition of (sparse) equilibrium systems.
subroutine cholsl(g, x, n)
Solution after decomposition.
subroutine chdec2(g, n, nrank, evmax, evmin, mon)
Cholesky decomposition (LARGE pos.
subroutine sort2k(a, n)
Quick sort 2.
subroutine devsig(n, diag, u, b, coef)
Calculate significances.
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
subroutine lltfwd(n, c, india, x)
Forward solution.
subroutine dbsvx(v, a, b, n)
Product symmetric matrix, vector.
subroutine equslvs(n, m, b, nm, c, india, l, x)
Solution of (sparse) equilibrium systems (after decomposition).
subroutine precons(p, n, b, nm, c, cu, a, l, s, nrkd)
Constrained (sparse) preconditioner, decomposition.
subroutine sort2i(a, n)
Quick sort 2 with index.
integer, parameter mpd
double precision
integer, parameter mpl
long integer
integer, parameter mps
single precision
integer, parameter mpi
integer
subroutine peend(icode, cmessage)
Print exit code.