 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-04-18 Time: 19:56:08 00004 00213 00214 00215 ! (1) Symmetric matrix W: decomposition, solution, inverse 00216 00225 00226 SUBROUTINE dchdec(w,n, aux) 00227 implicit none 00228 integer :: i 00229 integer :: ii 00230 integer :: j 00231 integer :: jj 00232 integer :: k 00233 integer :: kk 00234 integer :: l 00235 integer :: m 00236 00237 00238 DOUBLE PRECISION, INTENT(IN OUT) :: w(*) 00239 INTEGER, INTENT(IN) :: n 00240 DOUBLE PRECISION, INTENT(OUT) :: aux(n) 00241 00242 DOUBLE PRECISION :: b(*),x(*),v(*),suma,ratio 00243 ! ... 00244 DO i=1,n 00245 aux(i)=16.0D0*w((i*i+i)/2) ! save diagonal elements 00246 END DO 00247 ii=0 00248 DO i=1,n 00249 ii=ii+i 00250 IF(w(ii)+aux(i) /= aux(i)) THEN ! GT 00251 w(ii)=1.0D0/w(ii) ! (I,I) div ! 00252 ELSE 00253 w(ii)=0.0D0 00254 END IF 00255 jj=ii 00256 DO j=i+1,n 00257 ratio=w(i+jj)*w(ii) ! (I,J) (I,I) 00258 kk=jj 00259 DO k=j,n 00260 w(kk+j)=w(kk+j)-w(kk+i)*ratio ! (K,J) (K,I) 00261 kk=kk+k 00262 END DO ! K 00263 w(i+jj)=ratio ! (I,J) 00264 jj=jj+j 00265 END DO ! J 00266 END DO ! I 00267 00268 00269 RETURN 00270 00271 ENTRY dchslv(w,n,b, x) ! solution B -> X 00272 WRITE(*,*) 'before copy' 00273 DO i=1,n 00274 x(i)=b(i) 00275 END DO 00276 WRITE(*,*) 'after copy' 00277 ii=0 00278 DO i=1,n 00279 suma=x(i) 00280 DO k=1,i-1 00281 suma=suma-w(k+ii)*x(k) ! (K,I) 00282 END DO 00283 x(i)=suma 00284 ii=ii+i 00285 END DO 00286 WRITE(*,*) 'after forward' 00287 DO i=n,1,-1 00288 suma=x(i)*w(ii) ! (I,I) 00289 kk=ii 00290 DO k=i+1,n 00291 suma=suma-w(kk+i)*x(k) ! (K,I) 00292 kk=kk+k 00293 END DO 00294 x(i)=suma 00295 ii=ii-i 00296 END DO 00297 WRITE(*,*) 'after backward' 00298 RETURN 00299 00300 ENTRY dchinv(w,n,v) ! inversion 00301 ii=(n*n-n)/2 00302 DO i=n,1,-1 00303 suma=w(ii+i) ! (I,I) 00304 DO j=i,1,-1 00305 DO k=j+1,n 00306 l=MIN(i,k) 00307 m=MAX(i,k) 00308 suma=suma-w(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K) 00309 END DO 00310 v(ii+j)=suma ! (I,J) 00311 suma=0.0D0 00312 END DO 00313 ii=ii-i+1 00314 END DO 00315 END SUBROUTINE dchdec 00316 00323 00324 REAL FUNCTION condes(w,n,aux) 00325 implicit none 00326 real :: cond 00327 integer :: i 00328 integer :: ir 00329 integer :: is 00330 integer :: k 00331 integer :: kk 00332 real :: xla1 00333 real :: xlan 00334 00335 00336 DOUBLE PRECISION, INTENT(IN) :: w(*) 00337 INTEGER, INTENT(IN) :: n 00338 DOUBLE PRECISION, INTENT(IN OUT) :: aux(n) 00339 00340 DOUBLE PRECISION :: suma,sumu,sums 00341 00342 ir=1 00343 is=1 00344 DO i=1,n 00345 IF(w((i*i+i)/2) < w((is*is+is)/2)) is=i ! largest Dii 00346 IF(w((i*i+i)/2) > w((ir*ir+ir)/2)) ir=i ! smallest Dii 00347 END DO 00348 00349 sumu=0.0 ! find smallest eigenvalue 00350 sums=0.0 00351 DO i=n,1,-1 ! backward 00352 suma=0.0 00353 IF(i == ir) suma=1.0D0 00354 kk=(i*i+i)/2 00355 DO k=i+1,n 00356 suma=suma-w(kk+i)*aux(k) ! (K,I) 00357 kk=kk+k 00358 END DO 00359 aux(i)=suma 00360 sumu=sumu+aux(i)*aux(i) 00361 END DO 00362 xlan=SNGL(w((ir*ir+ir)/2)*DSQRT(sumu)) 00363 IF(xlan /= 0.0) xlan=1.0/xlan 00364 00365 DO i=1,n 00366 IF(i == is) THEN 00367 sums=1.0D0 00368 ELSE IF(i > is) THEN 00369 sums=sums+w(is+(i*i-i)/2)**2 00370 END IF 00371 END DO ! is Ws 00372 xla1=0.0 00373 IF(w((is*is+is)/2) /= 0.0) xla1=SNGL(DSQRT(sums)/w((is*is+is)/2)) 00374 00375 cond=0.0 00376 IF(xla1 > 0.0.AND.xlan > 0.0) cond=xla1/xlan 00377 ! estimated condition 00378 condes=cond 00379 END FUNCTION condes 00380 00381 00382 ! (2) Symmetric band matrix: decomposition, solution, complete 00383 ! inverse and band part of the inverse 00395 00396 SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M 00397 implicit none 00398 integer :: i 00399 integer :: j 00400 integer :: k 00401 ! M=MP1-1 N*M(M-1) dot operations 00402 00403 DOUBLE PRECISION, INTENT(IN OUT) :: w(mp1,n) 00404 INTEGER, INTENT(IN) :: mp1 00405 INTEGER, INTENT(IN) :: n 00406 DOUBLE PRECISION, INTENT(OUT) :: aux(n) 00407 ! decompos 00408 DOUBLE PRECISION :: v(mp1,n),b(n),x(n), vs(*),rxw 00409 ! ... 00410 DO i=1,n 00411 aux(i)=16.0D0*w(1,i) ! save diagonal elements 00412 END DO 00413 DO i=1,n 00414 IF(w(1,i)+aux(i) /= aux(i)) THEN 00415 w(1,i)=1.0/w(1,i) 00416 ELSE 00417 w(1,i)=0.0D0 ! singular 00418 END IF 00419 DO j=1,MIN(mp1-1,n-i) 00420 rxw=w(j+1,i)*w(1,i) 00421 DO k=1,MIN(mp1-1,n-i)+1-j 00422 w(k,i+j)=w(k,i+j)-w(k+j,i)*rxw 00423 END DO 00424 w(j+1,i)=rxw 00425 END DO 00426 END DO 00427 RETURN 00428 00429 ENTRY dbcslv(w,mp1,n,b, x) ! solution B -> X 00430 ! N*(2M-1) dot operations 00431 DO i=1,n 00432 x(i)=b(i) 00433 END DO 00434 DO i=1,n ! forward substitution 00435 DO j=1,MIN(mp1-1,n-i) 00436 x(j+i)=x(j+i)-w(j+1,i)*x(i) 00437 END DO 00438 END DO 00439 DO i=n,1,-1 ! backward substitution 00440 rxw=x(i)*w(1,i) 00441 DO j=1,MIN(mp1-1,n-i) 00442 rxw=rxw-w(j+1,i)*x(j+i) 00443 END DO 00444 x(i)=rxw 00445 END DO 00446 RETURN 00447 00448 ENTRY dbciel(w,mp1,n, v) ! V = inverse band matrix elements 00449 ! N*M*(M-1) dot operations 00450 DO i=n,1,-1 00451 rxw=w(1,i) 00452 DO j=i,MAX(1,i-mp1+1),-1 00453 DO k=j+1,MIN(n,j+mp1-1) 00454 rxw=rxw-v(1+ABS(i-k),MIN(i,k))*w(1+k-j,j) 00455 END DO 00456 v(1+i-j,j)=rxw 00457 rxw=0.0 00458 END DO 00459 END DO 00460 RETURN 00461 00462 ENTRY dbcinb(w,mp1,n, vs) ! VS = band part of inverse symmetric matrix 00463 ! N*M*(M-1) dot operations 00464 DO i=n,1,-1 00465 rxw=w(1,i) 00466 DO j=i,MAX(1,i-mp1+1),-1 00467 DO k=j+1,MIN(n,j+mp1-1) 00468 rxw=rxw-vs((MAX(i,k)*(MAX(i,k)-1))/2+MIN(i,k))*w(1+k-j,j) 00469 END DO 00470 vs((i*i-i)/2+j)=rxw 00471 rxw=0.0 00472 END DO 00473 ! DO J=MAX(1,I-MP1+1)-1,1,-1 00474 ! VS((I*I-I)/2+J)=0.0 00475 ! END DO 00476 END DO 00477 RETURN 00478 00479 ENTRY dbcinv(w,mp1,n, vs) ! V = inverse symmetric matrix 00480 ! N*N/2*(M-1) dot operations 00481 DO i=n,1,-1 00482 rxw=w(1,i) 00483 DO j=i,1,-1 00484 DO k=j+1,MIN(n,j+mp1-1) 00485 rxw=rxw-vs((MAX(i,k)*(MAX(i,k)-1))/2+MIN(i,k))*w(1+k-j,j) 00486 END DO 00487 vs((i*i-i)/2+j)=rxw 00488 rxw=0.0 00489 END DO 00490 END DO 00491 RETURN 00492 END SUBROUTINE dbcdec 00493 00500 00501 SUBROUTINE dbcprv(w,mp1,n,b) 00502 implicit none 00503 integer :: i 00504 integer :: j 00505 integer :: nj 00506 real :: rho 00507 00508 00509 DOUBLE PRECISION, INTENT(IN OUT) :: w(mp1,n) 00510 INTEGER, INTENT(IN) :: mp1 00511 INTEGER, INTENT(IN) :: n 00512 DOUBLE PRECISION, INTENT(OUT) :: b(n) 00513 00514 DOUBLE PRECISION :: ERR 00515 INTEGER :: irho(5) 00516 ! ... 00517 WRITE(*,*) ' ' 00518 WRITE(*,101) 00519 00520 DO i=1,n 00521 ERR=DSQRT(w(1,i)) 00522 nj=0 00523 DO j=2,mp1 00524 IF(i+1-j > 0.AND.nj < 5) THEN 00525 nj=nj+1 00526 rho=SNGL(w(j,i+1-j)/(ERR*DSQRT(w(1,i+1-j)))) 00527 irho(nj)=IFIX(100.0*ABS(rho)+0.5) 00528 IF(rho < 0.0) irho(nj)=-irho(nj) 00529 END IF 00530 END DO 00531 WRITE(*,102) i,b(i),ERR,(irho(j),j=1,nj) 00532 END DO 00533 WRITE(*,103) 00534 101 FORMAT(5X,'i Param',7X,'error',7X,' c(i,i-1) c(i,i-2)'/) 00535 102 FORMAT(1X,i5,2G12.4,1X,5I9) 00536 103 FORMAT(33X,'(correlation coefficients in percent)') 00537 END SUBROUTINE dbcprv 00538 00544 00545 SUBROUTINE dbcprb(w,mp1,n) 00546 implicit none 00547 integer :: i 00548 integer :: j 00549 00550 00551 DOUBLE PRECISION, INTENT(IN OUT) :: w(mp1,n) 00552 INTEGER, INTENT(IN) :: mp1 00553 INTEGER, INTENT(IN) :: n 00554 00555 00556 ! ... 00557 IF(mp1 > 6) RETURN 00558 WRITE(*,*) ' ' 00559 WRITE(*,101) 00560 DO i=1,n 00561 WRITE(*,102) i,(w(j,i),j=1,mp1) 00562 END DO 00563 WRITE(*,*) ' ' 00564 101 FORMAT(5X,'i Diag ') 00565 102 FORMAT(1X,i5,6G12.4) 00566 END SUBROUTINE dbcprb 00567 00568 00569 ! (3) Symmetric band matrix of band width m=1: decomposition, 00570 ! solution, complete, inverse and band part of the inverse 00571 00592 00593 SUBROUTINE db2dec(w,n, aux) 00594 implicit none 00595 integer :: i 00596 00597 00598 DOUBLE PRECISION, INTENT(IN OUT) :: w(2,n) 00599 INTEGER, INTENT(IN OUT) :: n 00600 DOUBLE PRECISION, INTENT(OUT) :: aux(n) 00601 00602 DOUBLE PRECISION :: v(2,n),b(n),x(n), rxw 00603 00604 DO i=1,n 00605 aux(i)=16.0D0*w(1,i) ! save diagonal elements 00606 END DO 00607 DO i=1,n-1 00608 IF(w(1,i)+aux(i) /= aux(i)) THEN 00609 w(1,i)=1.0D0/w(1,i) 00610 rxw=w(2,i)*w(1,i) 00611 w(1,i+1)=w(1,i+1)-w(2,i)*rxw 00612 w(2,i)=rxw 00613 ELSE ! singular 00614 w(1,i)=0.0D0 00615 w(2,i)=0.0D0 00616 END IF 00617 END DO 00618 IF(w(1,n)+aux(n) > aux(n)) THEN ! N 00619 w(1,n)=1.0D0/w(1,n) 00620 ELSE ! singular 00621 w(1,n)=0.0D0 00622 END IF 00623 RETURN 00624 00625 ENTRY db2slv(w,n,b, x) ! solution B -> X 00626 ! The equation W(original)*X=B is solved for X; input is B in vector X. 00627 DO i=1,n 00628 x(i)=b(i) 00629 END DO 00630 DO i=1,n-1 ! by forward substitution 00631 x(i+1)=x(i+1)-w(2,i)*x(i) 00632 END DO 00633 x(n)=x(n)*w(1,n) ! and backward substitution 00634 DO i=n-1,1,-1 00635 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1) 00636 END DO 00637 RETURN 00638 00639 ENTRY db2iel(w,n, v) ! V = inverse band matrix elements 00640 ! The band elements of the inverse of W(original) are calculated, 00641 ! and stored in V in the same order as in W. 00642 ! Remaining elements of the inverse are not calculated. 00643 v(1,n )= w(1,n) 00644 v(2,n-1)=-v(1,n )*w(2,n-1) 00645 DO i=n-1,3,-1 00646 v(1,i )= w(1,i )-v(2,i )*w(2,i ) 00647 v(2,i-1)=-v(1,i )*w(2,i-1) 00648 END DO 00649 v(1,2)= w(1,2)-v(2,2)*w(2,2) 00650 v(2,1)=-v(1,2)*w(2,1) 00651 v(1,1)= w(1,1)-v(2,1)*w(2,1) 00652 END SUBROUTINE db2dec 00653 00654 00655 ! (4) Symmetric band matrix of band width m=2: decomposition, 00656 ! solution, complete, inverse and band part of the inverse 00657 00678 00679 SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2) 00680 implicit none 00681 integer :: i 00682 00683 00684 DOUBLE PRECISION, INTENT(IN OUT) :: w(3,n) 00685 INTEGER, INTENT(IN OUT) :: n 00686 DOUBLE PRECISION, INTENT(OUT) :: aux(n) 00687 ! decompos 00688 00689 DOUBLE PRECISION :: v(3,n),b(n),x(n), rxw 00690 00691 DO i=1,n 00692 aux(i)=16.0D0*w(1,i) ! save diagonal elements 00693 END DO 00694 DO i=1,n-2 00695 IF(w(1,i)+aux(i) /= aux(i)) THEN 00696 w(1,i)=1.0D0/w(1,i) 00697 rxw=w(2,i)*w(1,i) 00698 w(1,i+1)=w(1,i+1)-w(2,i)*rxw 00699 w(2,i+1)=w(2,i+1)-w(3,i)*rxw 00700 w(2,i)=rxw 00701 rxw=w(3,i)*w(1,i) 00702 w(1,i+2)=w(1,i+2)-w(3,i)*rxw 00703 w(3,i)=rxw 00704 ELSE ! singular 00705 w(1,i)=0.0D0 00706 w(2,i)=0.0D0 00707 w(3,i)=0.0D0 00708 END IF 00709 END DO 00710 IF(w(1,n-1)+aux(n-1) > aux(n-1)) THEN 00711 w(1,n-1)=1.0D0/w(1,n-1) 00712 rxw=w(2,n-1)*w(1,n-1) 00713 w(1,n)=w(1,n)-w(2,n-1)*rxw 00714 w(2,n-1)=rxw 00715 ELSE ! singular 00716 w(1,n-1)=0.0D0 00717 w(2,n-1)=0.0D0 00718 END IF 00719 IF(w(1,n)+aux(n) > aux(n)) THEN 00720 w(1,n)=1.0D0/w(1,n) 00721 ELSE ! singular 00722 w(1,n)=0.0D0 00723 END IF 00724 RETURN 00725 00726 ENTRY db3slv(w,n,b, x) ! solution B -> X 00727 DO i=1,n 00728 x(i)=b(i) 00729 END DO 00730 DO i=1,n-2 ! by forward substitution 00731 x(i+1)=x(i+1)-w(2,i)*x(i) 00732 x(i+2)=x(i+2)-w(3,i)*x(i) 00733 END DO 00734 x(n)=x(n)-w(2,n-1)*x(n-1) 00735 x(n)=x(n)*w(1,n) ! and backward substitution 00736 x(n-1)=x(n-1)*w(1,n-1)-w(2,n-1)*x(n) 00737 DO i=n-2,1,-1 00738 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2) 00739 END DO 00740 RETURN 00741 00742 ENTRY db3iel(w,n, v) ! V = inverse band matrix elements 00743 ! The band elements of the inverse of W(original) are calculated, 00744 ! and stored in V in the same order as in W. 00745 ! Remaining elements of the inverse are not calculated. 00746 v(1,n )= w(1,n) 00747 v(2,n-1)=-v(1,n )*w(2,n-1) 00748 v(3,n-2)=-v(2,n-1)*w(2,n-2)-v(1,n )*w(3,n-2) 00749 v(1,n-1)= w(1,n-1) -v(2,n-1)*w(2,n-1) 00750 v(2,n-2)=-v(1,n-1)*w(2,n-2)-v(2,n-1)*w(3,n-2) 00751 v(3,n-3)=-v(2,n-2)*w(2,n-3)-v(1,n-1)*w(3,n-3) 00752 DO i=n-2,3,-1 00753 v(1,i )= w(1,i ) -v(2,i )*w(2,i )-v(3,i)*w(3,i ) 00754 v(2,i-1)=-v(1,i )*w(2,i-1)-v(2,i)*w(3,i-1) 00755 v(3,i-2)=-v(2,i-1)*w(2,i-2)-v(1,i)*w(3,i-2) 00756 END DO 00757 v(1,2)= w(1,2) -v(2,2)*w(2,2)-v(3,2)*w(3,2) 00758 v(2,1)=-v(1,2)*w(2,1)-v(2,2)*w(3,1) 00759 v(1,1)= w(1,1) -v(2,1)*w(2,1)-v(3,1)*w(3,1) 00760 END SUBROUTINE db3dec 00761 00762 00763 ! (5) Symmetric matrix with band structure, bordered by full row/col 00764 ! - is not yet included - 00765 00766 ! SUBROUTINE BSOLV1(N,CU,RU,CK,RK,CH, BK,UH, AU) ! 1 00767 ! Input: CU = 3*N array replaced by decomposition 00768 ! RU N array rhs 00769 ! CK diagonal element 00770 ! RK rhs 00771 ! CH N-vector 00772 00773 ! Aux: AU N-vector auxliliary array 00774 00775 ! Result: FK curvature 00776 ! BK variance 00777 ! UH smoothed data points 00778 00779 00780 ! DOUBLE PRECISION CU(3,N),CI(3,N),CK,BK,AU(N),UH(N) 00781 ! ... 00782 ! CALL BDADEC(CU,3,N, AU) ! decomposition 00783 ! CALL DBASLV(CU,3,N, RU,UH) ! solve for zero curvature 00784 ! CALL DBASLV(CU,3,N, CH,AU) ! solve for aux. vector 00785 ! CTZ=0.0D0 00786 ! ZRU=0.0D0 00787 ! DO I=1,N 00788 ! CTZ=CTZ+CH(I)*AU(I) ! cT z 00789 ! ZRU=ZRU+RY(I)*AU(I) ! zT ru 00790 ! END DO 00791 ! BK=1.0D0/(CK-CTZ) ! variance of curvature 00792 ! FK=BK *(RK-ZRU) ! curvature 00793 ! DO I=1,N 00794 ! UH(I)=UH(I)-FK*AU(I) ! smoothed data points 00795 ! END DO 00796 ! RETURN 00797 00798 ! ENTRY BINV1(N,CU,CI, FK,AU) 00799 ! DOUBLE PRECISION CI(3,N) 00800 ! ... 00801 ! CALL DBAIBM(CU,3,N, CI) ! block part of inverse 00802 ! DO I=1,N 00803 ! CI(1,I)=CI(1,I)+FK*AU(I)*AU(I) ! diagonal elements 00804 ! IF(I.LT.N) CI(2,I)=CI(2,I)+FK*AU(I)*AU(I+1) ! next diagonal 00805 ! IF(I.LT.N-1) CI(3,I)=CI(3,I)+FK*AU(I)*AU(I+2) ! next diagonal 00806 ! END DO 00807 00808 ! END 00809 00818 00819 SUBROUTINE dcfdec(w,n) 00820 00821 IMPLICIT NONE 00822 DOUBLE PRECISION, INTENT(OUT) :: w(*) 00823 INTEGER, INTENT(IN) :: n 00824 INTEGER :: i,j,k 00825 DOUBLE PRECISION :: epsm,gamm,xchi,beta,delta,theta 00826 00827 epsm=2.2D-16 ! machine precision 00828 gamm=0.0D0 ! max diagonal element 00829 xchi=0.0D0 ! max off-diagonal element 00830 DO k=1,n 00831 gamm=MAX(gamm,ABS(w((k*k+k)/2))) 00832 DO j=k+1,n 00833 xchi=MAX(xchi,ABS(w((j*j-j)/2+k))) 00834 END DO 00835 END DO 00836 beta=SQRT(MAX(gamm,xchi/MAX(1.0D0,SQRT(dfloat(n*n-1))),epsm)) 00837 delta=epsm*MAX(1.0D0,gamm+xchi) 00838 00839 DO k=1,n 00840 DO i=1,k-1 00841 w((k*k-k)/2+i)=w((k*k-k)/2+i)*w((i*i+i)/2) 00842 END DO 00843 DO j=k+1,n 00844 DO i=1,k-1 00845 w((j*j-j)/2+k)=w((j*j-j)/2+k)-w((k*k-k)/2+i)*w((j*j-j)/2+i) 00846 END DO 00847 END DO 00848 theta=0.0D0 00849 DO j=k+1,n 00850 theta=MAX(theta,ABS(w((j*j-j)/2+k))) 00851 END DO 00852 w((k*k+k)/2)=1.0D0/MAX(ABS(w((k*k+k)/2)),(theta/beta)**2,delta) 00853 DO j=k+1,n 00854 w((j*j+j)/2)=w((j*j+j)/2)-w((j*j-j)/2+k)**2*w((k*k+k)/2) 00855 END DO 00856 END DO ! K 00857 00858 END SUBROUTINE dcfdec 00859 00869 00870 SUBROUTINE dbfdec(w,mp1,n) 00871 00872 IMPLICIT NONE 00873 DOUBLE PRECISION, INTENT(OUT) :: w(mp1,n) 00874 INTEGER, INTENT(IN OUT) :: mp1 00875 INTEGER, INTENT(IN) :: n 00876 INTEGER :: i,j,k 00877 DOUBLE PRECISION :: epsm,gamm,xchi,beta,delta,theta 00878 00879 epsm=2.2D-16 ! machine precision 00880 gamm=0.0D0 ! max diagonal element 00881 xchi=0.0D0 ! max off-diagonal element 00882 DO k=1,n 00883 gamm=MAX(gamm,ABS(w(1,k))) 00884 DO j=2,MIN(mp1,n-k+1) 00885 xchi=MAX(xchi,ABS(w(j,k))) 00886 END DO 00887 END DO 00888 beta=SQRT(MAX(gamm,xchi/MAX(1.0D0,SQRT(dfloat(n*n-1))),epsm)) 00889 delta=epsm*MAX(1.0D0,gamm+xchi) 00890 00891 DO k=1,n 00892 DO i=2,MIN(mp1,k) 00893 w(i,k-i+1)=w(i,k-i+1)*w(1,k-i+1) 00894 END DO 00895 DO j=2,MIN(mp1,n-k+1) 00896 DO i=MAX(2,j+k+1-mp1),k 00897 w(j,k)=w(j,k)-w(k-i+2,i-1)*w(j-i+k+1,i-1) 00898 END DO 00899 END DO 00900 theta=0.0D0 00901 DO j=2,MIN(mp1,n-k+1) 00902 theta=MAX(theta,ABS(w(j,k))) 00903 END DO 00904 w(1,k)=1.0D0/MAX(ABS(w(1,k)),(theta/beta)**2,delta) 00905 DO j=2,MIN(mp1,n-k+1) 00906 w(1,k+j-1)=w(1,k+j-1)-w(1,k)*w(j,k)**2 00907 END DO 00908 END DO ! K 00909 00910 END SUBROUTINE dbfdec 00911 00912
1.8.0
