 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:07:48 00004 00051 00052 !---------------------------------------------------------------------- 00071 00072 SUBROUTINE sqminv(v,b,n,nrank,diag,next) ! matrix inversion 00073 IMPLICIT NONE 00074 INTEGER :: i 00075 INTEGER :: ij 00076 INTEGER :: j 00077 INTEGER :: jj 00078 INTEGER :: jk 00079 INTEGER :: jl 00080 INTEGER :: k 00081 INTEGER :: kk 00082 INTEGER :: l 00083 INTEGER :: last 00084 INTEGER :: lk 00085 INTEGER :: next0 00086 00087 DOUBLE PRECISION, INTENT(IN OUT) :: v(*) 00088 DOUBLE PRECISION, INTENT(OUT) :: b(n) 00089 INTEGER, INTENT(IN) :: n 00090 INTEGER, INTENT(OUT) :: nrank 00091 DOUBLE PRECISION, INTENT(OUT) :: diag(n) 00092 INTEGER, INTENT(OUT) :: next(n) 00093 DOUBLE PRECISION :: vkk 00094 DOUBLE PRECISION :: vjk 00095 00096 DOUBLE PRECISION, PARAMETER :: eps=1.0D-10 00097 ! ... 00098 next0=1 00099 l=1 00100 DO i=1,n 00101 next(i)=i+1 ! set "next" pointer 00102 diag(i)=ABS(v((i*i+i)/2)) ! save abs of diagonal elements 00103 END DO 00104 next(n)=-1 ! end flag 00105 00106 nrank=0 00107 DO i=1,n ! start of loop 00108 k =0 00109 vkk=0.0D0 00110 00111 j=next0 00112 last=0 00113 05 IF(j > 0) THEN 00114 jj=(j*j+j)/2 00115 IF(ABS(v(jj)) > MAX(ABS(vkk),eps*diag(j))) THEN 00116 vkk=v(jj) 00117 k=j 00118 l=last 00119 END IF 00120 last=j 00121 j=next(last) 00122 GO TO 05 00123 END IF 00124 00125 IF(k /= 0) THEN ! pivot found 00126 kk=(k*k+k)/2 00127 IF(l == 0) THEN 00128 next0=next(k) 00129 ELSE 00130 next(l)=next(k) 00131 END IF 00132 next(k)=0 ! index is used, reset 00133 nrank=nrank+1 ! increase rank and ... 00134 vkk =1.0/vkk 00135 v(kk) =-vkk 00136 b(k) =b(k)*vkk 00137 jk =kk-k 00138 jl =0 00139 DO j=1,n ! elimination 00140 IF(j == k) THEN 00141 jk=kk 00142 jl=jl+j 00143 ELSE 00144 IF(j < k) THEN 00145 jk=jk+1 00146 ELSE 00147 jk=jk+j-1 00148 END IF 00149 vjk =v(jk) 00150 v(jk)=vkk*vjk 00151 b(j) =b(j)-b(k)*vjk 00152 lk =kk-k 00153 DO l=1,j 00154 jl=jl+1 00155 IF(l == k) THEN 00156 lk=kk 00157 ELSE 00158 IF(l < k) THEN 00159 lk=lk+1 00160 ELSE 00161 lk=lk+l-1 00162 END IF 00163 v(jl)=v(jl)-v(lk)*vjk 00164 END IF 00165 END DO 00166 END IF 00167 END DO 00168 ELSE 00169 DO k=1,n 00170 IF(next(k) /= 0) THEN 00171 b(k)=0.0D0 ! clear vector element 00172 DO j=1,k 00173 IF(next(j) /= 0) v((k*k-k)/2+j)=0.0D0 ! clear matrix row/col 00174 END DO 00175 END IF 00176 END DO 00177 GO TO 10 00178 END IF 00179 END DO ! end of loop 00180 10 DO ij=1,(n*n+n)/2 00181 v(ij)=-v(ij) ! finally reverse sign of all matrix elements 00182 END DO 00183 END SUBROUTINE sqminv 00184 00197 00198 SUBROUTINE sqminl(v,b,n,nrank,diag,next) ! 00199 USE mpdef 00200 00201 IMPLICIT NONE 00202 INTEGER :: i 00203 INTEGER :: j 00204 INTEGER :: k 00205 INTEGER :: l 00206 INTEGER :: last 00207 INTEGER :: next0 00208 00209 DOUBLE PRECISION, INTENT(IN OUT) :: v(*) 00210 DOUBLE PRECISION, INTENT(OUT) :: b(n) 00211 INTEGER, INTENT(IN) :: n 00212 INTEGER, INTENT(OUT) :: nrank 00213 DOUBLE PRECISION, INTENT(OUT) :: diag(n) 00214 INTEGER, INTENT(OUT) :: next(n) 00215 INTEGER*8 :: i8 00216 INTEGER*8 :: j8 00217 INTEGER*8 :: jj 00218 INTEGER*8 :: k8 00219 INTEGER*8 :: kk 00220 INTEGER*8 :: kkmk 00221 INTEGER*8 :: jk 00222 INTEGER*8 :: jl 00223 INTEGER*8 :: llk 00224 INTEGER*8 :: ljl 00225 00226 DOUBLE PRECISION :: vkk 00227 DOUBLE PRECISION :: vjk 00228 00229 DOUBLE PRECISION, PARAMETER :: eps=1.0D-10 00230 ! ... 00231 next0=1 00232 l=1 00233 DO i=1,n 00234 i8=int8(i) 00235 next(i)=i+1 ! set "next" pointer 00236 diag(i)=ABS(v((i8*i8+i8)/2)) ! save abs of diagonal elements 00237 END DO 00238 next(n)=-1 ! end flag 00239 00240 nrank=0 00241 DO i=1,n ! start of loop 00242 k =0 00243 vkk=0.0D0 00244 j=next0 00245 last=0 00246 05 IF(j > 0) THEN 00247 j8=int8(j) 00248 jj=(j8*j8+j8)/2 00249 IF(ABS(v(jj)) > MAX(ABS(vkk),eps*diag(j))) THEN 00250 vkk=v(jj) 00251 k=j 00252 l=last 00253 END IF 00254 last=j 00255 j=next(last) 00256 GO TO 05 00257 END IF 00258 00259 IF(k /= 0) THEN ! pivot found 00260 k8=int8(k) 00261 kk=(k8*k8+k8)/2 00262 kkmk=kk-k8 00263 IF(l == 0) THEN 00264 next0=next(k) 00265 ELSE 00266 next(l)=next(k) 00267 END IF 00268 next(k)=0 ! index is used, reset 00269 nrank=nrank+1 ! increase rank and ... 00270 vkk =1.0/vkk 00271 v(kk) =-vkk 00272 b(k) =b(k)*vkk 00273 ! elimination 00274 jk =kkmk 00275 DO j=1,n 00276 IF(j == k) THEN 00277 jk=kk 00278 ELSE 00279 IF(j < k) THEN 00280 jk=jk+1 00281 ELSE 00282 jk=jk+int8(j)-1 00283 END IF 00284 v(jk)=v(jk)*vkk 00285 END IF 00286 END DO 00287 ! parallelize row loop 00288 ! slot of 128 'J' for next idle thread 00289 !$OMP PARALLEL DO & 00290 !$OMP PRIVATE(JL,JK,L,LJL,LLK,VJK,J8) & 00291 !$OMP SCHEDULE(DYNAMIC,128) 00292 DO j=n,1,-1 00293 j8=int8(j) 00294 jl=j8*(j8-1)/2 00295 IF(j /= k) THEN 00296 IF(j < k) THEN 00297 jk=kkmk+j8 00298 ELSE 00299 jk=k8+jl 00300 END IF 00301 vjk =v(jk)/vkk 00302 b(j) =b(j)-b(k)*vjk 00303 ljl=jl 00304 llk=kkmk 00305 DO l=1,MIN(j,k-1) 00306 ljl=ljl+1 00307 llk=llk+1 00308 v(ljl)=v(ljl)-v(llk)*vjk 00309 END DO 00310 ljl=ljl+1 00311 llk=kk 00312 DO l=k+1,j 00313 ljl=ljl+1 00314 llk=llk+l-1 00315 v(ljl)=v(ljl)-v(llk)*vjk 00316 END DO 00317 END IF 00318 END DO 00319 !$OMP END PARALLEL DO 00320 ELSE 00321 DO k=1,n 00322 k8=int8(k) 00323 kk=(k8*k8-k8)/2 00324 IF(next(k) /= 0) THEN 00325 b(k)=0.0D0 ! clear vector element 00326 DO j=1,k 00327 IF(next(j) /= 0) v(kk+int8(j))=0.0D0 ! clear matrix row/col 00328 END DO 00329 END IF 00330 END DO 00331 GO TO 10 00332 END IF 00333 END DO ! end of loop 00334 10 DO jj=1,(int8(n)*int8(n)+int8(n))/2 00335 v(jj)=-v(jj) ! finally reverse sign of all matrix elements 00336 END DO 00337 END SUBROUTINE sqminl 00338 00350 00351 SUBROUTINE devrot(n,diag,u,v,work,iwork) ! diagonalization 00352 IMPLICIT NONE 00353 00354 INTEGER, INTENT(IN) :: n 00355 DOUBLE PRECISION, INTENT(OUT) :: diag(n) 00356 DOUBLE PRECISION, INTENT(OUT) :: u(n,n) 00357 DOUBLE PRECISION, INTENT(IN) :: v(*) 00358 DOUBLE PRECISION, INTENT(OUT) :: work(n) 00359 INTEGER, INTENT(OUT) :: iwork(n) 00360 00361 00362 INTEGER, PARAMETER :: itmax=30 00363 DOUBLE PRECISION, PARAMETER :: tol=1.0D-16 00364 DOUBLE PRECISION, PARAMETER :: eps=1.0D-16 00365 00366 DOUBLE PRECISION :: f 00367 DOUBLE PRECISION :: g 00368 DOUBLE PRECISION :: h 00369 DOUBLE PRECISION :: sh 00370 DOUBLE PRECISION :: hh 00371 DOUBLE PRECISION :: b 00372 DOUBLE PRECISION :: p 00373 DOUBLE PRECISION :: r 00374 DOUBLE PRECISION :: s 00375 DOUBLE PRECISION :: c 00376 DOUBLE PRECISION :: workd 00377 00378 INTEGER :: ij 00379 INTEGER :: i 00380 INTEGER :: j 00381 INTEGER :: k 00382 INTEGER :: l 00383 INTEGER :: m 00384 INTEGER :: ll 00385 ! ... 00386 ! 1. part: symmetric matrix V reduced to tridiagonal from 00387 ij=0 00388 DO i=1,n 00389 DO j=1,i 00390 ij=ij+1 00391 u(i,j)=v(ij) ! copy half of symmetric matirx 00392 END DO 00393 END DO 00394 00395 DO i=n,2,-1 00396 l=i-2 00397 f=u(i,i-1) 00398 g=0.0D0 00399 IF(l /= 0) THEN 00400 DO k=1,l 00401 IF(ABS(u(i,k)) > tol) g=g+u(i,k)*u(i,k) 00402 END DO 00403 h=g+f*f 00404 END IF 00405 IF(g < tol) THEN ! G too small 00406 work(i)=f ! skip transformation 00407 h =0.0D0 00408 ELSE 00409 l=l+1 00410 sh=SQRT(h) 00411 IF(f >= 0.0D0) sh=-sh 00412 g=sh 00413 work(i)=sh 00414 h=h-f*g 00415 u(i,i-1)=f-g 00416 f=0.0D0 00417 DO j=1,l 00418 u(j,i)=u(i,j)/h 00419 g=0.0D0 00420 ! form element of a u 00421 DO k=1,j 00422 IF(ABS(u(j,k)) > tol.AND.ABS(u(i,k)) > tol) THEN 00423 g=g+u(j,k)*u(i,k) 00424 END IF 00425 END DO 00426 DO k=j+1,l 00427 IF(ABS(u(k,j)) > tol.AND.ABS(u(i,k)) > tol) THEN 00428 g=g+u(k,j)*u(i,k) 00429 END IF 00430 END DO 00431 work(j)=g/h 00432 f=f+g*u(j,i) 00433 END DO 00434 ! form k 00435 hh=f/(h+h) 00436 ! form reduced a 00437 DO j=1,l 00438 f=u(i,j) 00439 work(j)=work(j)-hh*f 00440 g=work(j) 00441 DO k=1,j 00442 u(j,k)=u(j,k)-f*work(k)-g*u(i,k) 00443 END DO 00444 END DO 00445 END IF 00446 diag(i)=h 00447 END DO 00448 00449 diag(1)=0.0D0 00450 work(1)=0.0D0 00451 00452 ! accumulation of transformation matrices 00453 DO i=1,n 00454 IF(diag(i) /= 0.0) THEN 00455 DO j=1,i-1 00456 g=0.0D0 00457 DO k=1,i-1 00458 g=g+u(i,k)*u(k,j) 00459 END DO 00460 DO k=1,i-1 00461 u(k,j)=u(k,j)-g*u(k,i) 00462 END DO 00463 END DO 00464 END IF 00465 diag(i)=u(i,i) 00466 u(i,i)=1.0D0 00467 DO j=1,i-1 00468 u(i,j)=0.0D0 00469 u(j,i)=0.0D0 00470 END DO 00471 END DO 00472 00473 ! 2. part: diagonalization of tridiagonal matrix 00474 DO i=2,n 00475 work(i-1)=work(i) 00476 END DO 00477 work(n)=0.0D0 00478 b=0.0D0 00479 f=0.0D0 00480 00481 DO l=1,n 00482 j=0 00483 h=eps*(ABS(diag(l))+ABS(work(l))) 00484 IF(b < h) b=h 00485 DO m=l,n 00486 IF(ABS(work(m)) <= b) GO TO 10 ! look for small sub-diagonal element 00487 END DO 00488 m=l 00489 10 IF(m == l) GO TO 30 00490 ! next iteration 00491 20 IF(j == itmax) THEN 00492 WRITE(*,*) 'DEVROT: Iteration limit reached' 00493 STOP 00494 END IF 00495 j=j+1 00496 g=diag(l) 00497 p=(diag(l+1)-g)/(2.0D0*work(l)) 00498 r=SQRT(1.0D0+p*p) 00499 diag(l)=work(l) 00500 IF(p < 0.0D0) diag(l)=diag(l)/(p-r) 00501 IF(p >= 0.0D0) diag(l)=diag(l)/(p+r) 00502 h=g-diag(l) 00503 DO i=l+1,n 00504 diag(i)=diag(i)-h 00505 END DO 00506 f=f+h 00507 ! QL transformation 00508 p=diag(m) 00509 c=1.0D0 00510 s=0.0D0 00511 DO i=m-1,l,-1 ! reverse loop 00512 g=c*work(i) 00513 h=c*p 00514 IF(ABS(p) >= ABS(work(i))) THEN 00515 c=work(i)/p 00516 r=SQRT(1.0D0+c*c) 00517 work(i+1)=s*p*r 00518 s=c/r 00519 c=1.0D0/r 00520 ELSE 00521 c=p/work(i) 00522 r=SQRT(1.0D0+c*c) 00523 work(i+1)=s*work(i)*r 00524 s=1.0D0/r 00525 c=c/r 00526 END IF 00527 p=c*diag(i)-s*g 00528 diag(i+1)=h+s*(c*g+s*diag(i)) 00529 ! form vector 00530 DO k=1,n 00531 h=u(k,i+1) 00532 u(k,i+1)=s*u(k,i)+c*h 00533 u(k,i)=c*u(k,i)-s*h 00534 END DO 00535 END DO 00536 work(l)=s*p 00537 diag(l)=c*p 00538 IF(ABS(work(l)) > b) GO TO 20 ! next iteration 00539 30 diag(l)=diag(l)+f 00540 END DO 00541 DO i=1,n 00542 iwork(i)=i 00543 END DO 00544 00545 m=1 00546 40 m=1+3*m ! determine initial increment 00547 IF(m <= n) GO TO 40 00548 50 m=m/3 00549 DO j=1,n-m ! sort with increment M 00550 l=j 00551 60 IF(diag(iwork(l+m)) > diag(iwork(l))) THEN ! compare 00552 ll=iwork(l+m) ! exchange the two index values 00553 iwork(l+m)=iwork(l) 00554 iwork(l)=ll 00555 l=l-m 00556 IF(l > 0) GO TO 60 00557 END IF 00558 END DO 00559 IF(m > 1) GO TO 50 00560 00561 DO i=1,n 00562 IF(iwork(i) /= i) THEN 00563 ! move vector from position I to the work area 00564 workd=diag(i) 00565 DO l=1,n 00566 work(l)=u(l,i) 00567 END DO 00568 k=i 00569 70 j=k 00570 k=iwork(j) 00571 iwork(j)=j 00572 IF(k /= i) THEN 00573 ! move vector from position K to the (free) position J 00574 diag(j)=diag(k) 00575 DO l=1,n 00576 u(l,j)=u(l,k) 00577 END DO 00578 GO TO 70 00579 END IF 00580 ! move vector from the work area to position J 00581 diag(j)=workd 00582 DO l=1,n 00583 u(l,j)=work(l) 00584 END DO 00585 END IF 00586 END DO 00587 END SUBROUTINE devrot 00588 00590 SUBROUTINE devsig(n,diag,u,b,coef) 00591 IMPLICIT NONE 00592 00593 INTEGER, INTENT(IN) :: n 00594 DOUBLE PRECISION, INTENT(IN) :: diag(n) 00595 DOUBLE PRECISION, INTENT(IN) :: u(n,n) 00596 DOUBLE PRECISION, INTENT(IN) :: b(n) 00597 DOUBLE PRECISION, INTENT(OUT) :: coef(n) 00598 INTEGER :: i 00599 INTEGER :: j 00600 DOUBLE PRECISION :: dsum 00601 ! ... 00602 DO i=1,n 00603 coef(i)=0.0D0 00604 IF(diag(i) > 0.0D0) THEN 00605 dsum=0.0D0 00606 DO j=1,n 00607 dsum=dsum+u(j,i)*b(j) 00608 END DO 00609 coef(i)=ABS(dsum)/SQRT(diag(i)) 00610 END IF 00611 END DO 00612 END SUBROUTINE devsig 00613 00614 00625 00626 SUBROUTINE devsol(n,diag,u,b,x,work) 00627 IMPLICIT NONE 00628 00629 INTEGER, INTENT(IN) :: n 00630 DOUBLE PRECISION, INTENT(IN) :: diag(n) 00631 DOUBLE PRECISION, INTENT(IN) :: u(n,n) 00632 DOUBLE PRECISION, INTENT(IN) :: b(n) 00633 DOUBLE PRECISION, INTENT(OUT) :: x(n) 00634 DOUBLE PRECISION, INTENT(OUT) :: work(n) 00635 INTEGER :: i 00636 INTEGER :: j 00637 INTEGER :: jj 00638 DOUBLE PRECISION :: s 00639 ! ... 00640 DO j=1,n 00641 s=0.0D0 00642 work(j)=0.0D0 00643 IF(diag(j) /= 0.0D0) THEN 00644 DO i=1,n 00645 ! j-th eigenvector is U(.,J) 00646 s=s+u(i,j)*b(i) 00647 END DO 00648 work(j)=s/diag(j) 00649 END IF 00650 END DO 00651 00652 DO j=1,n 00653 s=0.0D0 00654 DO jj=1,n 00655 s=s+u(j,jj)*work(jj) 00656 END DO 00657 x(j)=s 00658 END DO 00659 ! WRITE(*,*) 'DEVSOL' 00660 ! WRITE(*,*) 'X ',X 00661 END SUBROUTINE devsol 00662 00670 00671 SUBROUTINE devinv(n,diag,u,v) 00672 00673 IMPLICIT NONE 00674 INTEGER :: i 00675 INTEGER :: ij 00676 INTEGER :: j 00677 INTEGER :: k 00678 00679 INTEGER, INTENT(IN) :: n 00680 DOUBLE PRECISION, INTENT(IN) :: diag(n) 00681 DOUBLE PRECISION, INTENT(IN) :: u(n,n) 00682 DOUBLE PRECISION, INTENT(OUT) :: v(*) 00683 DOUBLE PRECISION:: dsum 00684 ! ... 00685 ij=0 00686 DO i=1,n 00687 DO j=1,i 00688 ij=ij+1 00689 dsum=0.0D0 00690 DO k=1,n 00691 IF(diag(k) /= 0.0D0) THEN 00692 dsum=dsum+u(i,k)*u(j,k)/diag(k) 00693 END IF 00694 END DO 00695 v(ij)=dsum 00696 END DO 00697 END DO 00698 END SUBROUTINE devinv 00699 00700 00719 SUBROUTINE choldc(g,n) 00720 IMPLICIT NONE 00721 INTEGER :: i 00722 INTEGER :: ii 00723 INTEGER :: j 00724 INTEGER :: jj 00725 INTEGER :: k 00726 INTEGER :: kk 00727 00728 DOUBLE PRECISION, INTENT(IN OUT) :: g(*) 00729 INTEGER, INTENT(IN) :: n 00730 00731 DOUBLE PRECISION :: ratio 00732 ! ... 00733 ii=0 00734 DO i=1,n 00735 ii=ii+i 00736 IF(g(ii) /= 0.0) g(ii)=1.0/g(ii) ! (I,I) div ! 00737 jj=ii 00738 DO j=i+1,n 00739 ratio=g(i+jj)*g(ii) ! (I,J) (I,I) 00740 kk=jj 00741 DO k=j,n 00742 g(kk+j)=g(kk+j)-g(kk+i)*ratio ! (K,J) (K,I) 00743 kk=kk+k 00744 END DO ! K 00745 g(i+jj)=ratio ! (I,J) 00746 jj=jj+j 00747 END DO ! J 00748 END DO ! I 00749 RETURN 00750 END SUBROUTINE choldc 00751 00763 SUBROUTINE cholsl(g,x,n) 00764 IMPLICIT NONE 00765 DOUBLE PRECISION :: dsum 00766 INTEGER :: i 00767 INTEGER :: ii 00768 INTEGER :: k 00769 INTEGER :: kk 00770 00771 DOUBLE PRECISION, INTENT(IN) :: g(*) 00772 DOUBLE PRECISION, INTENT(IN OUT) :: x(n) 00773 INTEGER, INTENT(IN) :: n 00774 00775 ii=0 00776 DO i=1,n 00777 dsum=x(i) 00778 DO k=1,i-1 00779 dsum=dsum-g(k+ii)*x(k) ! (K,I) 00780 END DO 00781 x(i)=dsum 00782 ii=ii+i 00783 END DO 00784 DO i=n,1,-1 00785 dsum=x(i)*g(ii) ! (I,I) 00786 kk=ii 00787 DO k=i+1,n 00788 dsum=dsum-g(kk+i)*x(k) ! (K,I) 00789 kk=kk+k 00790 END DO 00791 x(i)=dsum 00792 ii=ii-i 00793 END DO 00794 RETURN 00795 END SUBROUTINE cholsl 00796 00807 SUBROUTINE cholin(g,v,n) 00808 IMPLICIT NONE 00809 DOUBLE PRECISION :: dsum 00810 INTEGER :: i 00811 INTEGER :: ii 00812 INTEGER :: j 00813 INTEGER :: k 00814 INTEGER :: l 00815 INTEGER :: m 00816 00817 DOUBLE PRECISION, INTENT(IN) :: g(*) 00818 DOUBLE PRECISION, INTENT( OUT) :: v(*) 00819 INTEGER, INTENT(IN) :: n 00820 00821 ii=(n*n-n)/2 00822 DO i=n,1,-1 00823 dsum=g(ii+i) ! (I,I) 00824 DO j=i,1,-1 00825 DO k=j+1,n 00826 l=MIN(i,k) 00827 m=MAX(i,k) 00828 dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K) 00829 END DO 00830 v(ii+j)=dsum ! (I,J) 00831 dsum=0.0D0 00832 END DO 00833 ii=ii-i+1 00834 END DO 00835 END SUBROUTINE cholin 00836 00837 ! variable band matrix operations ---------------------------------- 00838 00859 00860 SUBROUTINE vabdec(n,val,ilptr) 00861 IMPLICIT NONE 00862 00863 INTEGER :: i 00864 INTEGER :: in 00865 INTEGER :: j 00866 INTEGER :: k 00867 INTEGER :: kj 00868 INTEGER :: mj 00869 INTEGER :: mk 00870 DOUBLE PRECISION :: sn 00871 DOUBLE PRECISION :: beta 00872 DOUBLE PRECISION :: delta 00873 DOUBLE PRECISION :: theta 00874 00875 INTEGER, INTENT(IN) :: n 00876 DOUBLE PRECISION, INTENT(IN OUT) :: val(*) 00877 INTEGER, INTENT(IN) :: ilptr(n) 00878 00879 DOUBLE PRECISION :: dgamma 00880 DOUBLE PRECISION :: xi 00881 DOUBLE PRECISION :: eps 00882 DOUBLE PRECISION :: prd 00883 DOUBLE PRECISION :: valkj 00884 00885 DOUBLE PRECISION, PARAMETER :: one=1.0D0 00886 DOUBLE PRECISION, PARAMETER :: two=2.0D0 00887 ! DATA EPS/0.0D0/ 00888 DATA eps/2.22044605E-16/ 00889 SAVE eps 00890 ! ... 00891 IF(eps == 0.0D0) THEN 00892 eps = two**(-12) 00893 10 eps = eps/two 00894 prd=one 00895 prd=prd+one*eps ! CALL dbaxpy(1,one,eps,prd) 00896 IF(prd > one) GO TO 10 00897 eps=eps*two ! EPS is machine presision 00898 WRITE(*,*) 'Machine precision is ',eps 00899 END IF 00900 00901 WRITE(*,*) 'Variable band matrix Cholesky decomposition' 00902 00903 dgamma=0.0D0 00904 i=1 00905 DO j=1,ilptr(n) ! loop thrugh all matrix elements 00906 IF(ilptr(i) == j) THEN ! diagonal element 00907 IF(val(j) <= 0.0D0) GO TO 01 ! exit loop for negative diag 00908 dgamma=MAX(dgamma,ABS(val(j))) ! max diagonal element 00909 i=i+1 00910 END IF 00911 END DO 00912 i=n+1 00913 01 in=i-1 ! IN positive diagonal elements 00914 WRITE(*,*) ' ',in,' positive diagonal elements' 00915 xi=0.0D0 00916 i=1 00917 DO j=1,ilptr(in) ! loop for positive diagonal elements 00918 ! through all matrix elements 00919 IF(ilptr(i) == j) THEN ! diagonal element 00920 i=i+1 00921 ELSE 00922 xi=MAX(xi,ABS(val(j))) ! Xi = abs(max) off-diagonal element 00923 END IF 00924 END DO 00925 00926 delta=eps*MAX(1.0D0,dgamma+xi) 00927 sn=1.0D0 00928 IF(n > 1) sn=1.0D0/SQRT(dfloat(n*n-1)) 00929 beta=SQRT(MAX(eps,dgamma,xi*sn)) ! beta 00930 WRITE(*,*) ' DELTA and BETA ',delta,beta 00931 00932 DO k=2,n 00933 mk=k-ilptr(k)+ilptr(k-1)+1 00934 00935 theta=0.0D0 00936 00937 DO j=mk,k 00938 mj=j-ilptr(j)+ilptr(j-1)+1 00939 kj=ilptr(k)-k+j ! index kj 00940 00941 DO i=MAX(mj,mk),j-1 00942 val(kj)=val(kj) & ! L_kj := L_kj - L_ki D_ii L_ji 00943 -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i) 00944 00945 END DO ! 00946 00947 theta=MAX(theta,ABS(val(kj))) ! maximum value of row 00948 00949 IF(j /= k) THEN 00950 IF(val(ilptr(j)) /= 0.0D0) THEN 00951 val(kj)=val(kj)/val(ilptr(j)) 00952 ELSE 00953 val(kj)=0.0D0 00954 END IF 00955 END IF ! L_kj := L_kj/D_jj ! D_kk 00956 00957 IF(j == k) THEN 00958 valkj=val(kj) 00959 IF(k <= in) THEN 00960 val(kj)=MAX(ABS(val(kj)),(theta/beta)**2,delta) 00961 IF(valkj /= val(kj)) THEN 00962 WRITE(*,*) ' Index K=',k 00963 WRITE(*,*) ' ',valkj,val(kj), (theta/beta)**2,delta,theta 00964 END IF 00965 END IF 00966 END IF 00967 END DO ! J 00968 00969 END DO ! K 00970 00971 DO k=1,n 00972 IF(val(ilptr(k)) /= 0.0D0) val(ilptr(k))=1.0D0/val(ilptr(k)) 00973 END DO 00974 00975 RETURN 00976 END SUBROUTINE vabdec 00977 00983 00984 SUBROUTINE vabmmm(n,val,ilptr) 00985 IMPLICIT NONE 00986 INTEGER :: k 00987 INTEGER :: kp 00988 INTEGER :: kr 00989 INTEGER :: ks 00990 INTEGER, INTENT(IN) :: n 00991 DOUBLE PRECISION, INTENT(IN OUT) :: val(*) 00992 INTEGER, INTENT(IN) :: ilptr(n) 00993 kr=1 00994 ks=1 00995 kp=1 00996 DO k=1,n 00997 IF(val(ilptr(k)) > val(ilptr(ks))) ks=k 00998 IF(val(ilptr(k)) < val(ilptr(kr))) kr=k 00999 IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k 01000 END DO 01001 WRITE(*,*) ' Index value ',ks,val(ilptr(ks)) 01002 WRITE(*,*) ' Index value ',kp,val(ilptr(kp)) 01003 WRITE(*,*) ' Index value ',kr,val(ilptr(kr)) 01004 01005 RETURN 01006 END SUBROUTINE vabmmm 01007 01018 01019 SUBROUTINE vabslv(n,val,ilptr,x) 01020 IMPLICIT NONE 01021 INTEGER :: j 01022 INTEGER :: k 01023 INTEGER :: mk 01024 01025 INTEGER, INTENT(IN) :: n 01026 DOUBLE PRECISION, INTENT(IN OUT) :: val(*) 01027 INTEGER, INTENT(IN) :: ilptr(n) 01028 DOUBLE PRECISION, INTENT(IN OUT) :: x(n) 01029 ! ... 01030 DO k=1,n ! forward loop 01031 mk=k-ilptr(k)+ilptr(k-1)+1 01032 DO j=mk,k-1 01033 x(k)=x(k)-val(ilptr(k)-k+j)*x(j) ! X_k := X_k - L_kj B_j 01034 END DO 01035 END DO ! K 01036 01037 DO k=1,n ! divide by diagonal elements 01038 x(k)=x(k)*val(ilptr(k)) ! X_k := X_k*D_kk 01039 END DO 01040 01041 DO k=n,1,-1 ! backward loop 01042 mk=k-ilptr(k)+ilptr(k-1)+1 01043 DO j=mk,k-1 01044 x(j)=x(j)-val(ilptr(k)-k+j)*x(k) ! X_j := X_j - L_kj X_k 01045 END DO 01046 END DO ! K 01047 END SUBROUTINE vabslv 01048 01049 ! matrix/vector products ------------------------------------------- 01050 01057 01058 DOUBLE PRECISION FUNCTION dbdot(n,dx,dy) 01059 IMPLICIT NONE 01060 INTEGER:: i 01061 DOUBLE PRECISION :: dtemp 01062 01063 INTEGER, INTENT(IN) :: n 01064 DOUBLE PRECISION, INTENT(IN) :: dx(*) 01065 DOUBLE PRECISION, INTENT(IN) :: dy(*) 01066 ! ... 01067 dtemp=0.0D0 01068 DO i = 1,MOD(n,5) 01069 dtemp=dtemp+dx(i)*dy(i) 01070 END DO 01071 DO i =MOD(n,5)+1,n,5 01072 dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) & 01073 +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4) 01074 END DO 01075 dbdot=dtemp 01076 END FUNCTION dbdot 01077 01086 01087 SUBROUTINE dbaxpy(n,da,dx,dy) 01088 IMPLICIT NONE 01089 INTEGER:: i 01090 01091 INTEGER, INTENT(IN) :: n 01092 DOUBLE PRECISION, INTENT(IN) :: dx(*) 01093 DOUBLE PRECISION, INTENT(IN OUT) :: dy(*) 01094 DOUBLE PRECISION, INTENT(IN) :: da 01095 ! ... 01096 DO i=1,MOD(n,4) 01097 dy(i)=dy(i)+da*dx(i) 01098 END DO 01099 DO i=MOD(n,4)+1,n,4 01100 dy(i )=dy(i )+da*dx(i ) 01101 dy(i+1)=dy(i+1)+da*dx(i+1) 01102 dy(i+2)=dy(i+2)+da*dx(i+2) 01103 dy(i+3)=dy(i+3)+da*dx(i+3) 01104 END DO 01105 END SUBROUTINE dbaxpy 01106 01115 01116 SUBROUTINE dbsvx(v,a,b,n) ! 01117 IMPLICIT NONE 01118 INTEGER:: i 01119 INTEGER:: ij 01120 INTEGER:: ijs 01121 INTEGER:: j 01122 01123 ! B := V * A 01124 ! N N*N N 01125 01126 INTEGER, INTENT(IN) :: n 01127 DOUBLE PRECISION, INTENT(IN) :: v(*) 01128 DOUBLE PRECISION, INTENT(IN) :: a(*) 01129 DOUBLE PRECISION, INTENT(OUT) :: b(*) 01130 01131 DOUBLE PRECISION:: dsum 01132 ijs=1 01133 DO i=1,n 01134 dsum=0.0 01135 ij=ijs 01136 DO j=1,n 01137 dsum=dsum+v(ij)*a(j) 01138 IF(j < i) THEN 01139 ij=ij+1 01140 ELSE 01141 ij=ij+j 01142 END IF 01143 END DO 01144 b(i)=dsum 01145 ijs=ijs+i 01146 END DO 01147 END SUBROUTINE dbsvx 01148 01157 01158 SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector 01159 IMPLICIT NONE 01160 INTEGER:: i 01161 INTEGER:: j 01162 01163 ! B := V * A 01164 ! N N*N N 01165 01166 INTEGER, INTENT(IN) :: n 01167 DOUBLE PRECISION, INTENT(IN) :: v(*) 01168 DOUBLE PRECISION, INTENT(IN) :: a(*) 01169 DOUBLE PRECISION, INTENT(OUT) :: b(*) 01170 01171 DOUBLE PRECISION:: dsum 01172 INTEGER*8 :: ij 01173 INTEGER*8 :: ijs 01174 ijs=1 01175 DO i=1,n 01176 dsum=0.0 01177 ij=ijs 01178 DO j=1,n 01179 dsum=dsum+v(ij)*a(j) 01180 IF(j < i) THEN 01181 ij=ij+1 01182 ELSE 01183 ij=ij+int8(j) 01184 END IF 01185 END DO 01186 b(i)=dsum 01187 ijs=ijs+int8(i) 01188 END DO 01189 END SUBROUTINE dbsvxl 01190 01198 01199 SUBROUTINE dbgax(a,x,y,m,n) 01200 IMPLICIT NONE 01201 INTEGER :: i 01202 INTEGER :: ij 01203 INTEGER :: j 01204 01205 DOUBLE PRECISION, INTENT(IN) :: a(*) 01206 DOUBLE PRECISION, INTENT(IN) :: x(*) 01207 DOUBLE PRECISION, INTENT(OUT) :: y(*) 01208 INTEGER, INTENT(IN) :: m 01209 INTEGER, INTENT(IN) :: n 01210 01211 ! ... 01212 ij=0 01213 DO i=1,m 01214 y(i)=0.0D0 01215 DO j=1,n 01216 ij=ij+1 01217 y(i)=y(i)+a(ij)*x(j) 01218 END DO 01219 END DO 01220 END SUBROUTINE dbgax 01221 01234 SUBROUTINE dbavat(v,a,w,n,ms) 01235 IMPLICIT NONE 01236 INTEGER :: i 01237 INTEGER :: ij 01238 INTEGER :: ijs 01239 INTEGER :: il 01240 INTEGER :: j 01241 INTEGER :: jk 01242 INTEGER :: k 01243 INTEGER :: l 01244 INTEGER :: lk 01245 INTEGER :: lkl 01246 INTEGER :: m 01247 01248 DOUBLE PRECISION, INTENT(IN) :: v(*) 01249 DOUBLE PRECISION, INTENT(IN) :: a(*) 01250 DOUBLE PRECISION, INTENT(INOUT) :: w(*) 01251 INTEGER, INTENT(IN) :: n 01252 INTEGER, INTENT(IN) :: ms 01253 01254 DOUBLE PRECISION :: cik 01255 ! ... 01256 m=ms 01257 IF (m > 0) THEN 01258 DO i=1,(m*m+m)/2 01259 w(i)=0.0D0 ! reset output matrix 01260 END DO 01261 ELSE 01262 m=-m 01263 END IF 01264 01265 il=-n 01266 ijs=0 01267 DO i=1,m ! do I 01268 ijs=ijs+i-1 ! 01269 il=il+n ! 01270 lkl=0 ! 01271 DO k=1,n ! do K 01272 cik=0.0D0 ! 01273 lkl=lkl+k-1 ! 01274 lk=lkl ! 01275 DO l=1,k ! do L 01276 lk=lk+1 ! . 01277 cik=cik+a(il+l)*v(lk) ! . 01278 END DO ! end do L 01279 DO l=k+1,n ! do L 01280 lk=lk+l-1 ! . 01281 cik=cik+a(il+l)*v(lk) ! . 01282 END DO ! end do L 01283 jk=k ! 01284 ij=ijs ! 01285 DO j=1,i ! do J 01286 ij=ij+1 ! . 01287 w(ij)=w(ij)+cik*a(jk) ! . 01288 jk=jk+n ! . 01289 END DO ! end do J 01290 END DO ! end do K 01291 END DO ! end do I 01292 END SUBROUTINE dbavat 01293 01303 01304 SUBROUTINE dbmprv(lun,x,v,n) 01305 IMPLICIT NONE 01306 INTEGER :: i 01307 INTEGER :: ii 01308 INTEGER :: ij 01309 INTEGER :: j 01310 INTEGER :: jj 01311 INTEGER :: l 01312 INTEGER :: m 01313 INTEGER :: mc(15) 01314 REAL :: pd 01315 REAL :: rho 01316 REAL :: err 01317 01318 INTEGER, INTENT(IN) :: lun 01319 DOUBLE PRECISION, INTENT(IN) :: x(*) 01320 DOUBLE PRECISION, INTENT(IN) :: v(*) 01321 INTEGER, INTENT(IN) :: n 01322 01323 WRITE(lun,103) 01324 WRITE(lun,101) 01325 ii=0 01326 DO i=1,n 01327 ij=ii 01328 ii=ii+i 01329 ERR=0.0 01330 IF(v(ii) > 0.0) ERR=SQRT(REAL(v(ii))) 01331 l=0 01332 jj=0 01333 DO j=1,i 01334 jj=jj+j 01335 ij=ij+1 01336 rho=0.0 01337 pd=REAL(v(ii)*v(jj)) 01338 IF(pd > 0.0) rho=REAL(v(ij))/SQRT(pd) 01339 l=l+1 01340 mc(l)=INT(100.0*ABS(rho)+0.5) 01341 IF(rho < 0.0) mc(l)=-mc(l) 01342 IF(j == i.OR.l == 15) THEN 01343 IF(j <= 15) THEN 01344 IF(j == i) THEN 01345 WRITE(lun,102) i,x(i),ERR,(mc(m),m=1,l-1) 01346 ELSE 01347 WRITE(lun,102) i,x(i),ERR,(mc(m),m=1,l) 01348 END IF 01349 ELSE 01350 IF(j == i) THEN 01351 WRITE(lun,103) (mc(m),m=1,l-1) 01352 ELSE 01353 WRITE(lun,103) (mc(m),m=1,l) 01354 END IF 01355 l=0 01356 END IF 01357 END IF 01358 END DO 01359 END DO 01360 WRITE(lun,104) 01361 ! 100 RETURN 01362 RETURN 01363 101 FORMAT(9X,'Param',7X,'error',7X,'correlation coefficients'/) 01364 102 FORMAT(1X,i5,2G12.4,1X,15I5) 01365 103 FORMAT(31X,15I5) 01366 104 FORMAT(33X,'(correlation coefficients in percent)') 01367 END SUBROUTINE dbmprv 01368 01376 01377 SUBROUTINE dbprv(lun,v,n) 01378 IMPLICIT NONE 01379 INTEGER, PARAMETER :: istp=6 01380 INTEGER :: i 01381 INTEGER :: ip 01382 INTEGER :: ipe 01383 INTEGER :: ipn 01384 INTEGER :: ips 01385 INTEGER :: k 01386 01387 INTEGER, INTENT(IN) :: lun 01388 DOUBLE PRECISION, INTENT(IN) :: v(*) 01389 INTEGER, INTENT(IN) :: n 01390 01391 WRITE(lun,101) 01392 01393 DO i=1,n 01394 ips=(i*i-i)/2 01395 ipe=ips+i 01396 ip =ips 01397 100 CONTINUE 01398 ipn=ip+istp 01399 WRITE(lun,102), i, ip+1-ips, (v(k),k=ip+1,MIN(ipn,ipe)) 01400 IF (ipn < ipe) THEN 01401 ip=ipn 01402 GO TO 100 01403 END IF 01404 END DO 01405 RETURN 01406 101 FORMAT(1X,'--- DBPRV -----------------------------------') 01407 102 FORMAT(1X,2I3,6G12.4) 01408 END SUBROUTINE dbprv 01409 01410 ! sort ------------------------------------------------------------- 01411 01418 01419 SUBROUTINE heapf(a,n) 01420 IMPLICIT NONE 01421 ! 01422 INTEGER :: i 01423 INTEGER :: j 01424 INTEGER :: l 01425 INTEGER :: r 01426 REAL :: at ! pivot key value 01427 01428 REAL, INTENT(IN OUT) :: a(*) 01429 INTEGER, INTENT(IN) :: n 01430 ! ... 01431 IF(n <= 1) RETURN 01432 l=n/2+1 01433 r=n 01434 10 IF(l > 1) THEN 01435 l=l-1 01436 at =a(l) 01437 ELSE 01438 at =a(r) 01439 a(r)=a(1) 01440 r=r-1 01441 IF(r == 1) THEN 01442 a(1)=at 01443 RETURN 01444 END IF 01445 END IF 01446 i=l 01447 j=l+l 01448 20 IF(j <= r) THEN 01449 IF(j < r) THEN 01450 IF(a(j) < a(j+1)) j=j+1 01451 END IF 01452 IF(at < a(j)) THEN 01453 a(i)=a(j) 01454 i=j 01455 j=j+j 01456 ELSE 01457 j=r+1 01458 END IF 01459 GO TO 20 01460 END IF 01461 a(i)=at 01462 GO TO 10 01463 END SUBROUTINE heapf 01464 01471 01472 SUBROUTINE sort1k(a,n) 01473 IMPLICIT NONE 01474 INTEGER :: nlev ! stack size 01475 PARAMETER (nlev=2*32) ! ... for N = 2**32 = 4.3 10**9 01476 INTEGER :: i 01477 INTEGER :: j 01478 INTEGER :: l 01479 INTEGER :: r 01480 INTEGER :: lev 01481 INTEGER :: lr(nlev) 01482 INTEGER :: lrh 01483 INTEGER :: maxlev 01484 INTEGER :: a1 ! pivot key 01485 INTEGER :: at ! pivot key 01486 01487 INTEGER, INTENT(IN OUT) :: a(*) 01488 INTEGER, INTENT(IN) :: n 01489 ! ... 01490 IF (n <= 0) RETURN 01491 maxlev=0 01492 lev=0 01493 l=1 01494 r=n 01495 10 IF(r-l == 1) THEN ! sort two elements L and R 01496 IF(a(l) > a(r)) THEN 01497 at=a(l) ! exchange L <-> R 01498 a(l)=a(r) 01499 a(r)=at 01500 END IF 01501 r=l 01502 END IF 01503 IF(r == l) THEN 01504 IF(lev <= 0) THEN 01505 ! WRITE(*,*) 'SORT1K (quicksort): maxlevel used/available =', 01506 ! + MAXLEV,'/64' 01507 RETURN 01508 END IF 01509 lev=lev-2 01510 l=lr(lev+1) 01511 r=lr(lev+2) 01512 ELSE 01513 ! LRH=(L+R)/2 01514 lrh=(l/2)+(r/2) ! avoid bit overflow 01515 IF(MOD(l,2) == 1.AND.MOD(r,2) == 1) lrh=lrh+1 01516 a1=a(lrh) ! middle 01517 i=l-1 ! find limits [J,I] with [L,R] 01518 j=r+1 01519 20 i=i+1 01520 IF(a(i) < a1) GO TO 20 01521 30 j=j-1 01522 IF(a(j) > a1) GO TO 30 01523 IF(i <= j) THEN 01524 at=a(i) ! exchange I <-> J 01525 a(i)=a(j) 01526 a(j)=at 01527 GO TO 20 01528 END IF 01529 IF(lev+2 > nlev) STOP 'SORT1K (quicksort): stack overflow' 01530 IF(r-i < j-l) THEN 01531 lr(lev+1)=l 01532 lr(lev+2)=j 01533 l=i 01534 ELSE 01535 lr(lev+1)=i 01536 lr(lev+2)=r 01537 r=j 01538 END IF 01539 lev=lev+2 01540 maxlev=MAX(maxlev,lev) 01541 END IF 01542 GO TO 10 01543 END SUBROUTINE sort1k 01544 01551 01552 SUBROUTINE sort2k(a,n) 01553 IMPLICIT NONE 01554 INTEGER:: nlev ! stack size 01555 PARAMETER (nlev=2*32) ! ... for N = 2**32 = 4.3 10**9 01556 INTEGER:: i 01557 INTEGER::j 01558 INTEGER::l 01559 INTEGER::r 01560 INTEGER::lev 01561 INTEGER::lr(nlev) 01562 INTEGER::lrh 01563 INTEGER::maxlev 01564 INTEGER::a1 ! pivot key 01565 INTEGER::a2 ! pivot key 01566 INTEGER::at ! pivot key 01567 01568 INTEGER, INTENT(IN OUT) :: a(2,*) 01569 INTEGER, INTENT(IN) :: n 01570 ! ... 01571 maxlev=0 01572 lev=0 01573 l=1 01574 r=n 01575 10 IF(r-l == 1) THEN ! sort two elements L and R 01576 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN 01577 at=a(1,l) ! exchange L <-> R 01578 a(1,l)=a(1,r) 01579 a(1,r)=at 01580 at=a(2,l) 01581 a(2,l)=a(2,r) 01582 a(2,r)=at 01583 END IF 01584 r=l 01585 END IF 01586 IF(r == l) THEN 01587 IF(lev <= 0) THEN 01588 WRITE(*,*) 'SORT2K (quicksort): maxlevel used/available =', maxlev,'/64' 01589 RETURN 01590 END IF 01591 lev=lev-2 01592 l=lr(lev+1) 01593 r=lr(lev+2) 01594 ELSE 01595 ! LRH=(L+R)/2 01596 lrh=(l/2)+(r/2) ! avoid bit overflow 01597 IF(MOD(l,2) == 1.AND.MOD(r,2) == 1) lrh=lrh+1 01598 a1=a(1,lrh) ! middle 01599 a2=a(2,lrh) 01600 i=l-1 ! find limits [J,I] with [L,R] 01601 j=r+1 01602 20 i=i+1 01603 IF(a(1,i) < a1) GO TO 20 01604 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20 01605 30 j=j-1 01606 IF(a(1,j) > a1) GO TO 30 01607 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30 01608 IF(i <= j) THEN 01609 at=a(1,i) ! exchange I <-> J 01610 a(1,i)=a(1,j) 01611 a(1,j)=at 01612 at=a(2,i) 01613 a(2,i)=a(2,j) 01614 a(2,j)=at 01615 GO TO 20 01616 END IF 01617 IF(lev+2 > nlev) STOP 'SORT2K (quicksort): stack overflow' 01618 IF(r-i < j-l) THEN 01619 lr(lev+1)=l 01620 lr(lev+2)=j 01621 l=i 01622 ELSE 01623 lr(lev+1)=i 01624 lr(lev+2)=r 01625 r=j 01626 END IF 01627 lev=lev+2 01628 maxlev=MAX(maxlev,lev) 01629 END IF 01630 GO TO 10 01631 END SUBROUTINE sort2k 01632 01640 01641 REAL FUNCTION chindl(n,nd) 01642 IMPLICIT NONE 01643 INTEGER :: m 01644 INTEGER, INTENT(IN) :: n 01645 INTEGER, INTENT(IN) :: nd 01646 ! 01647 REAL:: sn(3) 01648 REAL::table(30,3) 01649 ! REAL PN(3) 01650 ! DATA PN/0.31731,0.0455002785,2.69985E-3/ ! probabilities 01651 DATA sn/0.47523,1.690140,2.782170/ 01652 DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, & 01653 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, & 01654 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, & 01655 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, & 01656 4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, & 01657 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, & 01658 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, & 01659 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, & 01660 9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, & 01661 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, & 01662 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, & 01663 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/ 01664 SAVE sn,table 01665 ! ... 01666 IF(nd < 1) THEN 01667 chindl=0.0 01668 ELSE 01669 m=MAX(1,MIN(n,3)) ! 1, 2 or 3 sigmas 01670 IF(nd <= 30) THEN 01671 chindl=table(nd,m) ! from table 01672 ELSE ! approximation for ND > 30 01673 chindl=(sn(m)+SQRT(FLOAT(nd+nd-1)))**2/FLOAT(nd+nd) 01674 END IF 01675 END IF 01676 END FUNCTION chindl 01677 01711 01712 SUBROUTINE lltdec(n,c,india,nrkd) 01713 IMPLICIT NONE 01714 INTEGER :: i 01715 INTEGER :: j 01716 INTEGER :: k 01717 INTEGER :: kj 01718 INTEGER :: mj 01719 INTEGER :: mk 01720 DOUBLE PRECISION::diag 01721 01722 INTEGER, INTENT(IN) :: n 01723 DOUBLE PRECISION, INTENT(IN OUT) :: c(*) 01724 INTEGER, INTENT(IN) :: india(n) 01725 INTEGER, INTENT(OUT) :: nrkd 01726 01727 ! ... 01728 nrkd=0 01729 diag=0.0D0 01730 IF(c(india(1)) > 0.0) THEN 01731 c(india(1))=1.0D0/SQRT(c(india(1))) ! square root 01732 ELSE 01733 c(india(1))=0.0D0 01734 nrkd=-1 01735 END IF 01736 01737 DO k=2,n 01738 mk=k-india(k)+india(k-1)+1 ! first index in row K 01739 DO j=mk,k ! loop over row K with index J 01740 mj=1 01741 IF(j > 1) mj=j-india(j)+india(j-1)+1 ! first index in row J 01742 kj=india(k)-k+j ! index kj 01743 diag=c(india(j)) ! j-th diagonal element 01744 01745 DO i=MAX(mj,mk),j-1 01746 ! L_kj = L_kj - L_ki *D_ii *L_ji 01747 c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i) 01748 END DO ! I 01749 01750 IF(j /= k) c(kj)=c(kj)*diag 01751 END DO ! J 01752 01753 IF(diag+c(india(k)) > diag) THEN ! test for linear dependence 01754 c(india(k))=1.0D0/SQRT(c(india(k))) ! square root 01755 ELSE 01756 DO j=mk,k ! reset row K 01757 c(india(k)-k+j)=0.0D0 01758 END DO ! J 01759 IF(nrkd == 0) THEN 01760 nrkd=-k 01761 ELSE 01762 IF(nrkd < 0) nrkd=1 01763 nrkd=nrkd+1 01764 END IF 01765 END IF 01766 01767 END DO ! K 01768 RETURN 01769 END SUBROUTINE lltdec 01770 01771 01783 01784 SUBROUTINE lltfwd(n,c,india,x) 01785 IMPLICIT NONE 01786 INTEGER :: j 01787 INTEGER :: k 01788 01789 INTEGER, INTENT(IN) :: n 01790 DOUBLE PRECISION, INTENT(IN) :: c(*) 01791 INTEGER, INTENT(IN) :: india(n) 01792 DOUBLE PRECISION, INTENT(IN OUT) :: x(n) 01793 01794 x(1)=x(1)*c(india(1)) 01795 DO k=2,n ! forward loop 01796 DO j=k-india(k)+india(k-1)+1,k-1 01797 x(k)=x(k)-c(india(k)-k+j)*x(j) ! X_k := X_k - L_kj * B_j 01798 END DO ! J 01799 x(k)=x(k)*c(india(k)) 01800 END DO ! K 01801 RETURN 01802 END SUBROUTINE lltfwd 01803 01815 01816 SUBROUTINE lltbwd(n,c,india,x) 01817 IMPLICIT NONE 01818 INTEGER :: j 01819 INTEGER :: k 01820 01821 INTEGER, INTENT(IN) :: n 01822 DOUBLE PRECISION, INTENT(IN) :: c(*) 01823 INTEGER, INTENT(IN) :: india(n) 01824 DOUBLE PRECISION, INTENT(IN OUT) :: x(n) 01825 01826 DO k=n,2,-1 ! backward loop 01827 x(k)=x(k)*c(india(k)) 01828 DO j=k-india(k)+india(k-1)+1,k-1 01829 x(j)=x(j)-c(india(k)-k+j)*x(k) ! X_j := X_j - L_kj * X_k 01830 END DO ! J 01831 END DO ! K 01832 x(1)=x(1)*c(india(1)) 01833 END SUBROUTINE lltbwd 01834 01851 SUBROUTINE equdec(n,m,c,india,nrkd,nrkd2) 01852 IMPLICIT NONE 01853 REAL:: eps 01854 INTEGER:: i 01855 INTEGER:: j 01856 INTEGER:: jk 01857 INTEGER:: k 01858 INTEGER:: ntotal 01859 01860 INTEGER, INTENT(IN) :: n 01861 INTEGER, INTENT(IN) :: m 01862 DOUBLE PRECISION, INTENT(IN OUT) :: c(*) 01863 INTEGER, INTENT(IN OUT) :: india(n+m) 01864 INTEGER, INTENT(OUT) :: nrkd 01865 INTEGER, INTENT(OUT) :: nrkd2 01866 01867 01868 PARAMETER (eps=2.22044605E-16) 01869 ! ... 01870 ntotal=n+n*m+(m*m+m)/2 01871 01872 CALL lltdec(n,c,india,nrkd) ! decomposition G G^T 01873 DO i=1,m 01874 CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K 01875 END DO 01876 01877 jk=india(n)+n*m 01878 DO j=1,m 01879 DO k=1,j 01880 jk=jk+1 01881 c(jk)=0.0D0 ! product K K^T 01882 DO i=1,n 01883 c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i) 01884 END DO 01885 END DO 01886 END DO 01887 01888 india(n+1)=1 01889 DO i=2,m 01890 india(n+i)=india(n+i-1)+MIN(i,m) ! pointer for K K^T 01891 END DO 01892 01893 CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2) ! decomp. H H^T 01894 01895 ntotal=n+n*m+(m*m+m)/2 01896 01897 RETURN 01898 END SUBROUTINE equdec 01899 01915 SUBROUTINE equslv(n,m,c,india,x) ! solution vector 01916 IMPLICIT NONE 01917 INTEGER :: i 01918 INTEGER :: j 01919 01920 INTEGER, INTENT(IN) :: n 01921 INTEGER, INTENT(IN) :: m 01922 DOUBLE PRECISION, INTENT(IN) :: c(*) 01923 INTEGER, INTENT(IN) :: india(n+m) 01924 DOUBLE PRECISION, INTENT(IN OUT) :: x(n+m) 01925 01926 CALL lltfwd(n,c,india,x) ! result is u 01927 DO i=1,m 01928 DO j=1,n 01929 x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u 01930 END DO 01931 END DO 01932 CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v 01933 01934 01935 CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y 01936 DO i=1,m 01937 x(n+i)=-x(n+i) ! result is +y 01938 END DO 01939 01940 DO i=1,n 01941 DO j=1,m 01942 x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y 01943 END DO 01944 END DO 01945 CALL lltbwd(n,c,india,x) ! result is x 01946 END SUBROUTINE equslv 01947 01976 01977 SUBROUTINE precon(p,n,c,cu,a,s) 01978 IMPLICIT NONE 01979 INTEGER :: i 01980 INTEGER :: ii 01981 INTEGER :: j 01982 INTEGER :: jj 01983 INTEGER :: jk 01984 INTEGER :: k 01985 INTEGER :: kk 01986 01987 INTEGER, INTENT(IN) :: p 01988 INTEGER, INTENT(IN) :: n 01989 DOUBLE PRECISION, INTENT(IN) :: c(n) 01990 DOUBLE PRECISION, INTENT(OUT) :: cu(n) 01991 DOUBLE PRECISION, INTENT(IN OUT) :: a(n,p) 01992 DOUBLE PRECISION, INTENT(OUT) :: s((p*p+p)/2) 01993 01994 DOUBLE PRECISION :: div 01995 DOUBLE PRECISION :: ratio 01996 01997 DO i=1,(p*p+p)/2 01998 s(i)=0.0D0 01999 END DO 02000 DO i=1,n 02001 jk=0 02002 div=c(i) ! copy 02003 IF (div > 0.0D0) THEN 02004 cu(i)=1.0D0/DSQRT(div) 02005 ELSE 02006 cu(i)=0.0D0 02007 END IF 02008 DO j=1,p 02009 a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2} 02010 DO k=1,j 02011 jk=jk+1 02012 s(jk)=s(jk)+a(i,j)*a(i,k) ! S = symmetric matrix K K^T 02013 END DO 02014 END DO 02015 END DO 02016 02017 ii=0 02018 DO i=1,p ! S -> H D H^T (Cholesky) 02019 ii=ii+i 02020 IF(s(ii) /= 0.0D0) s(ii)=1.0D0/s(ii) 02021 jj=ii 02022 DO j=i+1,p 02023 ratio=s(i+jj)*s(ii) 02024 kk=jj 02025 DO k=j,p 02026 s(kk+j)=s(kk+j)-s(kk+i)*ratio 02027 kk=kk+k 02028 END DO ! K 02029 s(i+jj)=ratio 02030 jj=jj+j 02031 END DO ! J 02032 END DO ! I 02033 RETURN 02034 END SUBROUTINE precon 02035 02045 02046 SUBROUTINE presol(p,n,cu,a,s,x,y) ! solution 02047 IMPLICIT NONE 02048 INTEGER :: i 02049 INTEGER :: j 02050 INTEGER :: jj 02051 INTEGER :: k 02052 INTEGER :: kk 02053 02054 INTEGER, INTENT(IN) :: p 02055 INTEGER, INTENT(IN) :: n 02056 02057 DOUBLE PRECISION, INTENT(IN) :: cu(n) 02058 DOUBLE PRECISION, INTENT(IN) :: a(n,p) 02059 DOUBLE PRECISION, INTENT(IN) :: s((p*p+p)/2) 02060 DOUBLE PRECISION, INTENT(OUT) :: x(n+p) 02061 DOUBLE PRECISION, INTENT(IN) :: y(n+p) 02062 02063 DOUBLE PRECISION :: dsum 02064 02065 DO i=1,n+p 02066 x(i)=y(i) 02067 END DO 02068 DO i=1,n 02069 x(i)=x(i)*cu(i) ! u =C^{-1/2} y 02070 DO j=1,p 02071 x(n+j)=x(n+j)-a(i,j)*x(i) ! d - K u 02072 END DO 02073 END DO 02074 02075 jj=0 02076 DO j=1,p ! Cholesky solution for v 02077 dsum=x(n+j) 02078 DO k=1,j-1 02079 dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u 02080 END DO 02081 x(n+j)=dsum ! -> v 02082 jj=jj+j 02083 END DO 02084 02085 DO j=p,1,-1 ! solution for lambda 02086 dsum=x(n+j)*s(jj) 02087 kk=jj 02088 DO k=j+1,p 02089 dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v 02090 kk=kk+k 02091 END DO 02092 x(n+j)=-dsum ! -> lambda 02093 jj=jj-j 02094 END DO 02095 02096 DO i=1,n ! u - K^T lambda 02097 DO j=1,p 02098 x(i)=x(i)-a(i,j)*x(n+j) 02099 END DO 02100 END DO 02101 DO i=1,n 02102 x(i)=x(i)*cu(i) ! x = C^{-1/2} u 02103 END DO 02104 02105 END SUBROUTINE presol 02106 02107 02108 ! 090817 C. Kleinwort, DESY-FH1 02153 SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag) 02154 ! Double precision scratch arrays: 02155 ! VBND(N*(NBND+1)) = storage of band part 02156 ! VBDR(N* NBDR) = storage of border part 02157 ! AUX (N* NBDR) = intermediate results 02158 02159 ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only) 02160 02161 IMPLICIT NONE 02162 INTEGER:: i 02163 INTEGER:: ib 02164 INTEGER:: ij 02165 INTEGER:: ioff 02166 INTEGER:: ip 02167 INTEGER:: ip1 02168 INTEGER:: ip2 02169 INTEGER:: is 02170 INTEGER:: j 02171 INTEGER:: j0 02172 INTEGER:: jb 02173 INTEGER:: joff 02174 INTEGER:: mp1 02175 INTEGER:: nb1 02176 INTEGER:: nmb 02177 INTEGER:: npri 02178 INTEGER:: nrankb 02179 02180 DOUBLE PRECISION, INTENT(IN OUT) :: v(*) 02181 DOUBLE PRECISION, INTENT(OUT) :: b(n) 02182 INTEGER, INTENT(IN) :: n 02183 INTEGER, INTENT(IN) :: nbdr 02184 INTEGER, INTENT(IN) :: nbnd 02185 INTEGER, INTENT(IN) :: inv 02186 INTEGER, INTENT(OUT) :: nrank 02187 02188 DOUBLE PRECISION, INTENT(OUT) :: vbnd(n*(nbnd+1)) 02189 DOUBLE PRECISION, INTENT(OUT) :: vbdr(n*nbdr) 02190 DOUBLE PRECISION, INTENT(OUT) :: aux(n*nbdr) 02191 DOUBLE PRECISION, INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2) 02192 DOUBLE PRECISION, INTENT(OUT) :: vzru(nbdr) 02193 DOUBLE PRECISION, INTENT(OUT) :: scdiag(nbdr) 02194 INTEGER, INTENT(OUT) :: scflag(nbdr) 02195 02196 SAVE npri 02197 DATA npri / 100 / 02198 ! ... 02199 nrank=0 02200 nb1=nbdr+1 02201 mp1=nbnd+1 02202 nmb=n-nbdr 02203 ! copy band part 02204 DO i=nb1,n 02205 ip=(i*(i+1))/2 02206 is=0 02207 DO j=i,MIN(n,i+nbnd) 02208 ip=ip+is 02209 is=j 02210 ib=j-i+1 02211 vbnd(ib+(i-nb1)*mp1)=v(ip) 02212 END DO 02213 END DO 02214 ! copy border part 02215 IF (nbdr > 0) THEN 02216 ioff=0 02217 DO i=1,nbdr 02218 ip=(i*(i+1))/2 02219 is=0 02220 DO j=i,n 02221 ip=ip+is 02222 is=j 02223 vbdr(ioff+j)=v(ip) 02224 END DO 02225 ioff=ioff+n 02226 END DO 02227 END IF 02228 02229 CALL dbcdec(vbnd,mp1,nmb,aux) 02230 ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable 02231 ! CALL DBCPRB(VBND,MP1,NMB) 02232 ip=1 02233 DO i=1, nmb 02234 IF (vbnd(ip) <= 0.0D0) THEN 02235 npri=npri-1 02236 IF (npri >= 0) THEN 02237 IF (vbnd(ip) == 0.0D0) THEN 02238 PRINT *, ' SQMIBB matrix singular', n, nbdr, nbnd 02239 ELSE 02240 PRINT *, ' SQMIBB matrix not positive definite', n, nbdr, nbnd 02241 END IF 02242 END IF 02243 ! return zeros 02244 DO ip=1,n 02245 b(ip)=0.0D0 02246 END DO 02247 DO ip=1,(n*n+n)/2 02248 v(ip)=0.0D0 02249 END DO 02250 RETURN 02251 END IF 02252 ip=ip+mp1 02253 END DO 02254 nrank=nmb 02255 02256 IF (nbdr == 0) THEN ! special case NBDR=0 02257 02258 CALL dbcslv(vbnd,mp1,nmb,b,b) 02259 IF (inv > 0) THEN 02260 IF (inv > 1) THEN 02261 CALL dbcinv(vbnd,mp1,nmb,v) 02262 ELSE 02263 CALL dbcinb(vbnd,mp1,nmb,v) 02264 END IF 02265 END IF 02266 02267 ELSE ! general case NBDR>0 02268 02269 ioff=nb1 02270 DO ib=1,nbdr 02271 ! solve for aux. vectors 02272 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff)) 02273 ! zT ru 02274 vzru(ib)=b(ib) 02275 DO i=0,nmb-1 02276 vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i) 02277 END DO 02278 ioff=ioff+n 02279 END DO 02280 ! solve for band part only 02281 CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1)) 02282 ! Ck - cT z 02283 ip=0 02284 ioff=nb1 02285 DO ib=1,nbdr 02286 joff=nb1 02287 DO jb=1,ib 02288 ip=ip+1 02289 vbk(ip)=v(ip) 02290 DO i=0,nmb-1 02291 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i) 02292 END DO 02293 joff=joff+n 02294 END DO 02295 ioff=ioff+n 02296 END DO 02297 ! solve border part 02298 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag) 02299 IF (nrankb == nbdr) THEN 02300 nrank=nrank+nbdr 02301 ELSE 02302 npri=npri-1 02303 IF (npri >= 0) PRINT *, ' SQMIBB undef border ', n, nbdr, nbnd, nrankb 02304 DO ib=1,nbdr 02305 vzru(ib)=0.0D0 02306 END DO 02307 DO ip=(nbdr*nbdr+nbdr)/2,1,-1 02308 vbk(ip)=0.0D0 02309 END DO 02310 END IF 02311 ! smoothed data points 02312 ioff=nb1 02313 DO ib=1, nbdr 02314 b(ib) = vzru(ib) 02315 DO i=0,nmb-1 02316 b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i) 02317 END DO 02318 ioff=ioff+n 02319 END DO 02320 ! inverse requested ? 02321 IF (inv > 0) THEN 02322 IF (inv > 1) THEN 02323 CALL dbcinv(vbnd,mp1,nmb,v) 02324 ELSE 02325 CALL dbcinb(vbnd,mp1,nmb,v) 02326 END IF 02327 ! expand/correct from NMB to N 02328 ip1=(nmb*nmb+nmb)/2 02329 ip2=(n*n+n)/2 02330 DO i=nmb-1,0,-1 02331 j0=0 02332 IF (inv == 1) j0=MAX(0,i-nbnd) 02333 DO j=i,j0,-1 02334 v(ip2)=v(ip1) 02335 ioff=nb1 02336 DO ib=1,nbdr 02337 joff=nb1 02338 DO jb=1,nbdr 02339 ij=MAX(ib,jb) 02340 ij=(ij*ij-ij)/2+MIN(ib,jb) 02341 v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j) 02342 joff=joff+n 02343 END DO 02344 ioff=ioff+n 02345 END DO 02346 ip1=ip1-1 02347 ip2=ip2-1 02348 END DO 02349 ip1=ip1-j0 02350 ip2=ip2-j0 02351 02352 DO ib=nbdr,1,-1 02353 v(ip2)=0.0D0 02354 joff=nb1 02355 DO jb=1,nbdr 02356 ij=MAX(ib,jb) 02357 ij=(ij*ij-ij)/2+MIN(ib,jb) 02358 v(ip2)=v(ip2)-vbk(ij)*aux(i+joff) 02359 joff=joff+n 02360 END DO 02361 ip2=ip2-1 02362 END DO 02363 END DO 02364 02365 DO ip=(nbdr*nbdr+nbdr)/2,1,-1 02366 v(ip2)=vbk(ip) 02367 ip2=ip2-1 02368 END DO 02369 02370 END IF 02371 END IF 02372 02373 END SUBROUTINE sqmibb
1.8.0
