|
N.I.M.R.O.D. |  |
|
N.I.M.R.O.D. | |
00001 00002 ! ------------------------------------------------ 00003 ! Cholesky decomposition. 00004 00005 ! input n size of matrix 00006 ! input A Symmetric positive def. matrix 00007 ! output aa lower deomposed matrix 00008 ! uses choldc1(int,MAT,VEC) 00009 ! ------------------------------------------------ 00010 Subroutine choldc(n,A,aa) 00011 integer n 00012 real*8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1) 00013 integer i,j, ialloc 00014 real*8, pointer :: p(:) 00015 allocate(p(0:n-1),stat=ialloc) 00016 00017 aa = A 00018 00019 call choldc1(n, aa, p) 00020 00021 do i = 0, n-1 00022 aa(i,i) = p(i) 00023 do j = i + 1, n-1 00024 aa(i,j) = 0.d0 00025 end do 00026 end do 00027 deallocate(p) 00028 return 00029 End 00030 00031 ! ----------------------------------------------------- 00032 ! Inverse of Cholesky decomposition. 00033 00034 ! input n size of matrix 00035 ! input A Symmetric positive def. matrix 00036 ! output aa inverse of lower decomposed matrix 00037 ! uses choldc1(int,MAT,VEC) 00038 ! ----------------------------------------------------- 00039 Subroutine choldcsl(n,A,aa) 00040 integer n 00041 real*8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1) 00042 integer i,j,k, ialloc 00043 real*8 sum 00044 real*8, pointer :: p(:) 00045 allocate(p(0:n-1),stat=ialloc) 00046 00047 aa = A 00048 00049 call choldc1(n, aa, p) 00050 00051 do i = 0, n-1 00052 aa(i,i) = 1.d0 / p(i) 00053 do j = i + 1, n-1 00054 sum = 0.d0 00055 do k = i, j-1 00056 sum = sum - aa(j,k) * aa(k,i) 00057 end do 00058 aa(j,i) = sum / p(j) 00059 end do 00060 end do 00061 deallocate(p) 00062 return 00063 End 00064 00065 ! ---------------------------------------------------------------------- 00066 ! Computation of Determinant of the matrix using Cholesky decomposition 00067 00068 ! input n size of matrix 00069 ! input a Symmetric positive def. matrix 00070 ! return det(a) 00071 ! uses choldc(int,MAT,MAT) 00072 ! ---------------------------------------------------------------------- 00073 real*8 Function choldet(n,a) 00074 integer n 00075 real*8 a(0:n-1,0:n-1) 00076 real*8, pointer :: c(:,:) 00077 real*8 d 00078 integer i, ialloc 00079 allocate(c(0:n-1,0:n-1),stat=ialloc) 00080 d=1.d0 00081 call choldc(n,a,c) 00082 do i = 0, n-1 00083 d = d * c(i,i) 00084 end do 00085 choldet = d * d 00086 deallocate(c) 00087 return 00088 End 00089 00090 00091 ! --------------------------------------------------- 00092 ! Matrix inverse using Cholesky decomposition 00093 00094 ! input n size of matrix 00095 ! input A Symmetric positive def. matrix 00096 ! output aa inverse of A 00097 ! uses choldc1(int,MAT,VEC) 00098 ! --------------------------------------------------- 00099 Subroutine cholsl(n,A,aa) 00100 integer n 00101 real*8 A(0:n-1,0:n-1), aa(0:n-1,0:n-1) 00102 integer i,j,k 00103 00104 call choldcsl(n,A,aa) 00105 00106 do i = 0, n-1 00107 do j = i + 1, n-1 00108 aa(i,j) = 0.d0 00109 end do 00110 end do 00111 00112 do i = 0, n-1 00113 aa(i,i) = aa(i,i) * aa(i,i) 00114 do k = i + 1, n-1 00115 aa(i,i) = aa(i,i) + aa(k,i) * aa(k,i) 00116 end do 00117 00118 do j = i + 1, n-1 00119 do k = j, n-1 00120 aa(i,j) = aa(i,j) + aa(k,i) * aa(k,j) 00121 end do 00122 end do 00123 end do 00124 do i = 0, n-1 00125 do j = 0, i-1 00126 aa(i,j) = aa(j,i) 00127 end do 00128 end do 00129 return 00130 End 00131 00132 ! ------------------------------------------------- 00133 ! main method for Cholesky decomposition. 00134 ! 00135 ! input n size of matrix 00136 ! input/output a matrix 00137 ! output p vector of resulting diag of a 00138 ! author: <Vadum Kutsyy, kutsyy@hotmail.com> 00139 ! ------------------------------------------------- 00140 Subroutine choldc1(n,a,p) 00141 integer n 00142 real*8 a(0:n-1,0:n-1), p(0:n-1) 00143 integer i,j,k 00144 real*8 sum 00145 do i = 0, n-1 00146 do j = i, n-1 00147 sum = a(i,j) 00148 do k = i - 1, 0, -1 00149 sum = sum - a(i,k) * a(j,k) 00150 end do 00151 if (i.eq.j) then 00152 if (sum <= 0.d0) & 00153 print *,' the matrix is not positive definite!' 00154 p(i) = dsqrt(sum) 00155 else 00156 a(j,i) = sum / p(i) 00157 end if 00158 end do 00159 end do 00160 return 00161 End 00162 00163 ! print a square real matrix A of size n with caption s 00164 ! (n items per line). 00165 Subroutine MatPrint(s,n,A) 00166 character*(*) s 00167 integer n 00168 real*8 A(0:n-1,0:n-1) 00169 integer i,j 00170 print *,' ' 00171 print *,' ', s 00172 do i=0, n-1 00173 write(*,10) (A(i,j),j=0,n-1) 00174 end do 00175 return 00176 10 format(10F10.6) 00177 End 00178 00179 Integer Function Check_Matrix(n,A) 00180 integer n 00181 real*8 A(0:n-1,0:n-1) 00182 integer i,j,k 00183 real*8 sum 00184 Check_Matrix=1 00185 do i = 0, n-1 00186 do j = i, n-1 00187 sum = A(i,j) 00188 do k = i - 1, 0, -1 00189 sum = sum - a(i,k) * a(j,k) 00190 end do 00191 if (i.eq.j) then 00192 if (sum <= 0.d0) Check_Matrix=0 00193 end if 00194 end do 00195 end do 00196 return 00197 End 00198 00199 !******************************************* 00200 !* MULTIPLICATION OF TWO SQUARE REAL * 00201 !* MATRICES * 00202 !* --------------------------------------- * 00203 !* INPUTS: A MATRIX N*N * 00204 !* B MATRIX N*N * 00205 !* N INTEGER * 00206 !* --------------------------------------- * 00207 !* OUTPUTS: C MATRIX N*N PRODUCT A*B * 00208 !* * 00209 !******************************************* 00210 Subroutine MATMULT(n,A,B,C) 00211 integer n 00212 real*8 A(0:n-1,0:n-1), B(0:n-1,0:n-1), C(0:n-1,0:n-1) 00213 real*8 SUM 00214 integer I,J,K 00215 do I=0, n-1 00216 do J=0, n-1 00217 SUM= 0.d0 00218 do K=0, n-1 00219 SUM=SUM+A(I,K)*B(K,J) 00220 end do 00221 C(I,J)=SUM 00222 end do 00223 end do 00224 return 00225 End 00226 00227 00228 !------------------------------------------------------------------------------ 00229 ! 00230 ! VERSION : 1.0 00231 ! 00232 ! MODULE: tools.f90 function PHID 00255 00256 00257 DOUBLE PRECISION FUNCTION PHID(Z) 00258 !* 00259 !* Normal distribution probabilities accurate to 1D-15. 00260 !* Z = no. of standard deviations from the mean. 00261 !* 00262 !* Based upon algorithm 5666 for the error function, from: 00263 !* Hart, J.F. et al, 'Computer Approximations', Wiley 1968 00264 !* 00265 !* Programmer: Alan Miller 00266 !* 00267 !* Latest revision - 30 March 1986 00268 !* 00269 DOUBLE PRECISION, INTENT(IN)::Z 00270 DOUBLE PRECISION P0, P1, P2, P3, P4, P5, P6, 00271 Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7, 00272 P, EXPNTL, CUTOFF, ROOTPI, ZABS 00273 PARAMETER(& 00274 P0 = 220.2068679123761D0, & 00275 P1 = 221.2135961699311D0, & 00276 P2 = 112.0792914978709D0,& 00277 P3 = 33.91286607838300D0,& 00278 P4 = 6.373962203531650D0,& 00279 P5 = .7003830644436881D0, & 00280 P6 = .03526249659989109D0) 00281 PARAMETER(& 00282 Q0 = 440.4137358247522D0,& 00283 Q1 = 793.8265125199484D0, & 00284 Q2 = 637.3336333788311D0,& 00285 Q3 = 296.5642487796737D0, & 00286 Q4 = 86.78073220294608D0,& 00287 Q5 = 16.06417757920695D0, & 00288 Q6 = 1.755667163182642D0,& 00289 Q7 = .08838834764831844D0) 00290 PARAMETER( ROOTPI = 2.506628274631001D0 ) 00291 PARAMETER( CUTOFF = 7.071067811865475D0 ) 00292 !* 00293 ZABS = ABS(Z) 00294 !* 00295 !* |Z| > 37 00296 !* 00297 IF (ZABS .GT. 37) THEN 00298 P = 0 00299 ELSE 00300 !* 00301 !* |Z| <= 37 00302 !* 00303 EXPNTL = EXP(-ZABS**2/2) 00304 !* 00305 !* |Z| < CUTOFF = 10/SQRT(2) 00306 !* 00307 IF ( ZABS .LT. CUTOFF ) THEN 00308 P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS & 00309 + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS & 00310 + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS & 00311 + Q0) 00312 !* 00313 !* |Z| >= CUTOFF. 00314 !* 00315 ELSE 00316 P = EXPNTL/(ZABS + 1/(ZABS + 2/(ZABS + 3/(ZABS + 4/ & 00317 (ZABS + 0.65D0)))))/ROOTPI 00318 END IF 00319 END IF 00320 IF ( Z .GT. 0 ) P = 1 - P 00321 PHID = P 00322 END 00323 00324 !------------------------------------------------------------------------------ 00325 ! 00326 ! VERSION : 1.0 00327 ! 00328 ! MODULE: tools.f90 function maxt 00348 subroutine dmaxt(maxt,delta,m) 00349 00350 implicit none 00351 00352 integer,intent(in)::m 00353 double precision,dimension(m),intent(in)::delta 00354 double precision,intent(out)::maxt 00355 integer::i 00356 00357 maxt=Dabs(delta(1)) 00358 do i=2,m 00359 if(Dabs(delta(i)).gt.maxt)then 00360 maxt=Dabs(delta(i)) 00361 end if 00362 end do 00363 00364 return 00365 end subroutine dmaxt 00366 00367 !------------------------------------------------------------------------------ 00368 ! 00369 ! VERSION : 1.0 00370 ! 00371 ! MODULE: tools.f90 function dchole 00393 subroutine dchole(a,k,nq,idpos) 00394 00395 implicit none 00396 00397 integer,intent(in)::k,nq 00398 integer,intent(inout)::idpos 00399 double precision,dimension(k*(k+3)/2),intent(inout)::a 00400 00401 integer::i,ii,i1,i2,i3,m,j,k2,jmk 00402 integer::ijm,irm,jji,jjj,l,jj,iil,jjl,il 00403 integer,dimension(k)::is 00404 double precision ::term,xn,diag,p 00405 equivalence (term,xn) 00406 00407 00408 ! ss programme de resolution d'un systeme lineaire symetrique 00409 ! 00410 ! k ordre du systeme / 00411 ! nq nombre de seconds membres 00412 ! 00413 ! en sortie les seconds membres sont remplaces par les solutions 00414 ! correspondantes 00415 ! 00416 00417 i2=0 00418 ii=0 00419 idpos=0 00420 k2=k+nq 00421 ! calcul des elements de la matrice 00422 do i=1,k 00423 ii=i*(i+1)/2 00424 ! elements diagonaux 00425 diag=a(ii) 00426 i1=ii-i 00427 if(i-1.ne.0) goto 1 00428 if(i-1.eq.0) goto 4 00429 1 i2=i-1 00430 do l=1,i2 00431 m=i1+l 00432 p=a(m) 00433 p=p*p 00434 if(is(l).lt.0) goto 2 00435 if(is(l).ge.0) goto 3 00436 2 p=-p 00437 3 diag=diag-p 00438 end do 00439 00440 4 if(diag.lt.0) goto 5 00441 if(diag.eq.0) goto 50 00442 if(diag.gt.0) goto 6 00443 5 is(i)=-1 00444 idpos=idpos+1 00445 diag=-dsqrt(-diag) 00446 a(ii)=-diag 00447 goto 7 00448 6 is(i)=1 00449 diag=dsqrt(diag) 00450 a(ii)=diag 00451 ! elements non diagonaux 00452 7 i3=i+1 00453 do j=i3,k2 00454 jj=j*(j-1)/2+i 00455 jmk=j-k-1 00456 if(jmk.le.0) goto 9 00457 if(jmk.gt.0) goto 8 00458 8 jj=jj-jmk*(jmk+1)/2 00459 9 term=a(jj) 00460 if(i-1.ne.0) goto 10 00461 if(i-1.eq.0) goto 13 00462 10 do l=1,i2 00463 iil=ii-l 00464 jjl=jj-l 00465 p=a(iil)*a(jjl) 00466 il=i-l 00467 if(is(il).lt.0) goto 11 00468 if(is(il).ge.0) goto 12 00469 11 p=-p 00470 12 term=term-p 00471 end do 00472 13 a(jj)=term/diag 00473 end do 00474 end do 00475 00476 ! calcul des solutions 00477 jj=ii-k+1 00478 do l=1,nq 00479 jj=jj+k 00480 i=k-1 00481 14 jji=jj+i 00482 xn=a(jji) 00483 if(i-k+1.lt.0) goto 20 00484 if(i-k+1.ge.0) goto 22 00485 20 j=k-1 00486 21 jjj=jj+j 00487 ijm=i+1+j*(j+1)/2 00488 xn=xn-a(jjj)*a(ijm) 00489 if(j-i-1.le.0) goto 22 00490 if(j-i-1.gt.0) goto 30 00491 30 j=j-1 00492 goto 21 00493 22 irm=(i+1)*(i+2)/2 00494 a(jji)=xn/a(irm) 00495 if(i.le.0) cycle 00496 if(i.gt.0) goto 40 00497 40 i=i-1 00498 go to 14 00499 end do 00500 50 continue 00501 return 00502 end subroutine dchole 00503 00504 00505 !------------------------------------------------------------------------------ 00506 ! 00507 ! VERSION : 1.0 00508 ! 00509 ! MODULE: tools.f90 function dmfsd 00531 subroutine dmfsd(a,n,eps,ier) 00532 ! 00533 ! FACTORISATION DE CHOLESKY D'UNE MATRICE SDP 00534 ! MATRICE = TRANSPOSEE(T)*T 00535 ! ENTREE : TABLEAU A CONTENANT LA PARTIE SUPERIEURE STOCKEE COLONNE 00536 ! PAR COLONNE DE LA METRICE A FACTORISER 00537 ! SORTIE : A CONTIENT LA PARTIE SUPPERIEURE DE LA MATRICE triangulaire T 00538 ! 00539 ! SUBROUTINE APPELE PAR DSINV 00540 ! 00541 ! N : DIM. MATRICE 00542 ! EPS : SEUIL DE TOLERANCE 00543 ! IER = 0 PAS D'ERREUR 00544 ! IER = -1 ERREUR 00545 ! IER = K COMPRIS ENTRE 1 ET N, WARNING, LE CALCUL CONTINUE 00546 ! 00547 implicit none 00548 00549 integer,intent(in)::n 00550 integer,intent(inout)::ier 00551 double precision,intent(in)::eps 00552 double precision,dimension(n*(n+1)/2),intent(inout)::A 00553 double precision :: dpiv,dsum,tol 00554 integer::i,k,l,kpiv,ind,lend,lanf,lind 00555 00556 ! 00557 ! TEST ON WRONG INPUT PARAMETER N 00558 ! 00559 dpiv=0.d0 00560 if (n-1.lt.0) goto 12 00561 if (n-1.ge.0) ier=0 00562 ! 00563 ! INITIALIZE DIAGONAL-LOOP 00564 ! 00565 kpiv=0 00566 do k=1,n 00567 kpiv=kpiv+k 00568 ind=kpiv 00569 lend=k-1 00570 ! 00571 ! CALCULATE TOLERANCE 00572 ! 00573 tol=dabs(eps*sngl(A(kpiv))) 00574 ! 00575 ! START FACTORIZATION-LOOP OVER K-TH ROW 00576 ! 00577 do i=k,n 00578 dsum=0.d0 00579 if (lend.lt.0) goto 2 00580 if (lend.eq.0) goto 4 00581 if (lend.gt.0) goto 2 00582 ! 00583 ! START INNER LOOP 00584 ! 00585 2 do l=1,lend 00586 lanf=kpiv-l 00587 lind=ind-l 00588 dsum=dsum+A(lanf)*A(lind) 00589 end do 00590 00591 ! 00592 ! END OF INNEF LOOP 00593 ! 00594 ! TRANSFORM ELEMENT A(IND) 00595 ! 00596 4 dsum=A(ind)-dsum 00597 if (i-k.ne.0) goto 10 00598 if (i-k.eq.0) goto 5 00599 ! TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE 00600 ! 00601 00602 00603 5 if (sngl(dsum)-tol.le.0) goto 6 00604 if (sngl(dsum)-tol.gt.0) goto 9 00605 6 if (dsum.le.0) goto 12 00606 if (dsum.gt.0) goto 7 00607 7 if (ier.le.0) goto 8 00608 if (ier.gt.0) goto 9 00609 8 ier=k-1 00610 ! 00611 ! COMPUTE PIVOT ELEMENT 00612 ! 00613 9 dpiv=dsqrt(dsum) 00614 A(kpiv)=dpiv 00615 dpiv=1.D0/dpiv 00616 goto 11 00617 ! 00618 ! CALCULATE TERMS IN ROW 00619 ! 00620 10 A(ind)=dsum*dpiv 00621 11 ind=ind+i 00622 end do 00623 end do 00624 00625 ! 00626 ! END OF DIAGONAL-LOOP 00627 ! 00628 return 00629 12 ier=-1 00630 return 00631 00632 end subroutine dmfsd 00633 00634 00635 00636 !------------------------------------------------------------------------------ 00637 ! 00638 ! VERSION : 1.0 00639 ! 00640 ! MODULE: tools.f90 function dsinv 00663 subroutine dsinv(A,N,EPS,IER,DET) 00664 00665 ! 00666 ! INVERSION D'UNE MATRICE SYMETRIQUE DEFINIE POSITIVE : 00667 ! 00668 ! MATRICE = TRANSPOSEE(T)*T 00669 ! INERSE(MATRICE) = INVERSE(T)*INVERSE(TRANSPOSEE(T)) 00670 ! 00671 ! A : TABLEAU CONTENANT LA PARTIE SUPERIEURE DE LA MATRICE A INVERSER 00672 ! STOCKEE COLONNE PAR COLONNE 00673 ! DIM. MATRICE A INVERSER = N 00674 ! DIM. TABLEAU A = N*(N+1)/2 00675 ! 00676 ! EPS : SEUIL DE TOLERANCE AU-DESSOUS DUQUEL UN PIVOT EST CONSIDERE 00677 ! COMME NUL 00678 ! 00679 ! IER : CODE D'ERREUR 00680 ! IER=0 PAS D'ERREUR 00681 ! IER=-1 ERREUR SUR LA DIM.N OU MATRICE PAS DEFINIE POSITIVE 00682 ! IER=1 PERTE DE SIGNIFICANCE, LE CALCUL CONTINUE 00683 ! 00684 implicit none 00685 00686 integer,intent(in)::n 00687 integer,intent(inout)::ier 00688 double precision,intent(inout)::eps 00689 double precision,intent(inout),optional::det 00690 double precision,dimension(n*(n+1)/2),intent(inout)::A 00691 double precision::din,work 00692 integer::ind,ipiv,i,j,k,l,min,kend,lhor,lver,lanf 00693 00694 ! 00695 ! FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMFSD 00696 ! A=TRANSPOSE(T) * T 00697 ! 00698 00699 call dmfsd(A,n,eps,ier) 00700 00701 det=0.d0 00702 00703 if (ier.lt.0) goto 9 00704 if (ier.ge.0) det=0.d0 00705 ! 00706 ! INVERT UPPER TRIANGULAR MATRIX T 00707 ! PREPARE INVERSION-LOOP 00708 ! 00709 ! 00710 ! calcul du log du determinant 00711 00712 do i=1,n 00713 det=det+dlog(A(i*(i+1)/2)) 00714 end do 00715 det=2*det 00716 ipiv=n*(n+1)/2 00717 ind=ipiv 00718 ! 00719 ! INITIALIZE INVERSION-LOOP 00720 ! 00721 do i=1,n 00722 din=1.d0/A(ipiv) 00723 A(ipiv)=din 00724 min=n 00725 kend=i-1 00726 lanf=n-kend 00727 if (kend.le.0) goto 5 00728 if (kend.gt.0) j=ind 00729 ! 00730 ! INITIALIZE ROW-LOOP 00731 ! 00732 do k=1,kend 00733 work=0.d0 00734 min=min-1 00735 lhor=ipiv 00736 lver=j 00737 ! 00738 ! START INNER LOOP 00739 ! 00740 do l=lanf,min 00741 lver=lver+1 00742 lhor=lhor+l 00743 work=work+A(lver)*A(lhor) 00744 end do 00745 ! 00746 ! END OF INNER LOOP 00747 ! 00748 A(j)=-work*din 00749 j=j-min 00750 end do 00751 00752 ! 00753 ! END OF ROW-LOOP 00754 ! 00755 5 ipiv=ipiv-min 00756 ind=ind-1 00757 end do 00758 00759 ! 00760 ! END OF INVERSION-LOOP 00761 ! 00762 ! CALCULATE INVERSE(A) BY MEANS OF INVERSE(T) 00763 ! INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T)) 00764 ! INITIALIZE MULTIPLICATION-LOOP 00765 ! 00766 do i=1,n 00767 ipiv=ipiv+i 00768 j=ipiv 00769 ! 00770 ! INITIALIZE ROW-LOOP 00771 ! 00772 do k=i,n 00773 work=0.d0 00774 lhor=j 00775 ! 00776 ! START INNER LOOP 00777 ! 00778 do l=k,n 00779 lver=lhor+k-i 00780 work=work+A(lhor)*A(lver) 00781 lhor=lhor+l 00782 end do 00783 ! 00784 ! END OF INNER LOOP 00785 ! 00786 A(j)=work 00787 j=j+k 00788 end do 00789 end do 00790 00791 ! 00792 ! END OF ROW-AND MULTIPLICATION-LOOP 00793 ! 00794 9 return 00795 end subroutine dsinv 00796 00797 00798 !------------------------------------------------------------------------------ 00799 ! 00800 ! VERSION : 1.0 00801 ! 00802 ! MODULE: tools.f90 function gonfle 00822 subroutine inflateDiag(v,m) 00823 use ToleranceThreshold 00824 implicit none 00825 integer, intent(in)::m 00826 double precision,dimension(m*(m+3)/2),intent(inout) :: v 00827 00828 integer :: i,nql,ii,nfmax,idpos,ncount 00829 double precision,dimension(m*(m+3)/2) :: fu 00830 double precision :: da,dm,ga,tr 00831 00832 nql=0 00833 nfmax=m*(m+1)/2 00834 tr=0.d0 00835 da=gonfle_da 00836 dm=gonfle_dm 00837 do i=1,m 00838 ii=i*(i+1)/2 00839 tr=tr+dabs(v(ii)) 00840 end do 00841 tr=tr/dble(m) 00842 ncount=0 00843 ga=gonfle_ga 00844 400 continue 00845 do i=1,nfmax+m 00846 fu(i)=v(i) 00847 end do 00848 do i=1,m 00849 ii=i*(i+1)/2 00850 if (v(ii).ne.0) then 00851 fu(ii)=v(ii)+da*((1.d0-ga)*dabs(v(ii))+ga*tr) 00852 else 00853 fu(ii)=da*ga*tr 00854 endif 00855 end do 00856 call dchole(fu,m,nql,idpos) 00857 if (idpos.ne.0) then 00858 ncount=ncount+1 00859 if (ncount.le.3.or.ga.ge.1.d0) then 00860 da=da*dm 00861 else 00862 ga=ga*dm 00863 if (ga.gt.1.d0) ga=1.d0 00864 endif 00865 if (ncount > 100000) then 00866 fu=0.d0 00867 do i=1,m 00868 ii=i*(i+1)/2 00869 fu(ii)=1 00870 end do 00871 return 00872 end if 00873 goto 400 00874 else 00875 v=fu 00876 return 00877 end if 00878 end subroutine inflateDiag 00879 00880 00881 !------------------------------------------------------------------------------ 00882 ! 00883 ! VERSION : 1.0 00884 ! 00885 ! MODULE: tools.f90 function FINDInv 00907 SUBROUTINE FINDInv(matrix, inverse, n, errorflag) 00908 IMPLICIT NONE 00909 !Declarations 00910 INTEGER, INTENT(IN) :: n 00911 INTEGER, INTENT(OUT) :: errorflag !Return error status. -1 for error, 0 for normal 00912 double precision, INTENT(IN), DIMENSION(n,n) :: matrix !Input matrix 00913 double precision, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix 00914 00915 LOGICAL :: FLAG = .TRUE. 00916 INTEGER :: i, j, k, l 00917 double precision :: m 00918 double precision, DIMENSION(n,2*n) :: augmatrix !augmented matrix 00919 00920 !Augment input matrix with an identity matrix 00921 DO i = 1, n 00922 DO j = 1, 2*n 00923 IF (j <= n ) THEN 00924 augmatrix(i,j) = matrix(i,j) 00925 ELSE IF ((i+n) == j) THEN 00926 augmatrix(i,j) = 1 00927 Else 00928 augmatrix(i,j) = 0 00929 ENDIF 00930 END DO 00931 END DO 00932 00933 !Reduce augmented matrix to upper traingular form 00934 DO k =1, n-1 00935 IF (augmatrix(k,k) == 0) THEN 00936 FLAG = .FALSE. 00937 DO i = k+1, n 00938 IF (augmatrix(i,k) /= 0) THEN 00939 DO j = 1,2*n 00940 augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j) 00941 END DO 00942 FLAG = .TRUE. 00943 EXIT 00944 ENDIF 00945 IF (FLAG .EQV. .FALSE.) THEN 00946 PRINT*, "Matrix is non - invertible" 00947 inverse = 0 00948 errorflag = -1 00949 return 00950 ENDIF 00951 END DO 00952 ENDIF 00953 DO j = k+1, n 00954 m = augmatrix(j,k)/augmatrix(k,k) 00955 DO i = k, 2*n 00956 augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i) 00957 END DO 00958 END DO 00959 END DO 00960 00961 !Test for invertibility 00962 DO i = 1, n 00963 IF (augmatrix(i,i) == 0) THEN 00964 PRINT*, "Matrix is non - invertible" 00965 inverse = 0 00966 errorflag = -1 00967 return 00968 ENDIF 00969 END DO 00970 00971 !Make diagonal elements as 1 00972 DO i = 1 , n 00973 m = augmatrix(i,i) 00974 DO j = i , (2 * n) 00975 augmatrix(i,j) = (augmatrix(i,j) / m) 00976 END DO 00977 END DO 00978 00979 !Reduced right side half of augmented matrix to identity matrix 00980 DO k = n-1, 1, -1 00981 DO i =1, k 00982 m = augmatrix(i,k+1) 00983 DO j = k, (2*n) 00984 augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m 00985 END DO 00986 END DO 00987 END DO 00988 00989 !store answer 00990 DO i =1, n 00991 DO j = 1, n 00992 inverse(i,j) = augmatrix(i,j+n) 00993 END DO 00994 END DO 00995 errorflag = 0 00996 END SUBROUTINE FINDinv 00997 00998 !------------------------------------------------------------------------------ 00999 ! 01000 ! VERSION : 1.0 01001 ! 01002 ! MODULE: tools.f90 module m_vpropre 01031 module m_vpropre 01032 01033 contains 01034 SUBROUTINE EIGEN(NM,N,A,WR,WI,VECS,D,INTC) 01035 01036 !C 01037 !C CALCUL DES VALEURS ET DES VECTEURS PROPRES DE LA MATRICE DE TRANSITION 01038 !C 01039 INTEGER N,NM,LINF,LSUP,INTC(NM),i,j,k,l 01040 DOUBLE PRECISION A(NM,NM),WR(NM),WI(NM),VECS(NM,NM),D(NM) 01041 !C 01042 !C 01043 !C AMELIORATION DE LA PRECISION 01044 !C 01045 CALL BALANC(NM,N,A,LINF,LSUP,D) 01046 !C 01047 !C TRANSFORMATION DE LA MATRICE SOUS LA FORME 01048 !C HESSENBERG 01049 !C 01050 CALL ELMHES(NM,N,LINF,LSUP,A,INTC) 01051 CALL ELTRAN(NM,N,LINF,LSUP,A,INTC,VECS) 01052 !C 01053 !C CALCUL VALEURS PROPRES ET VECTEURS PROPRES 01054 !C VLR = PARTIE REELLE 01055 !C VLI = PARTIE IMAGINAIRE 01056 !C 01057 CALL HQR2(NM,N,LINF,LSUP,A,VECS,WR,WI) 01058 !C 01059 !C ALCUL DES VECTEURS PROPRES DE LA MATRICE ORIGINELLE 01060 !C 01061 CALL BALBAK(NM,N,LINF,LSUP,D,VECS) 01062 01063 RETURN 01064 END SUBROUTINE EIGEN 01065 01066 !******************************************************************** 01067 01068 01069 subroutine balanc(nm,n,a,linf,lsup,d) 01070 ! 01071 integer n,nm,linf,lsup 01072 double precision a(nm,nm),d(nm),radix,sqrdx 01073 01074 integer i,j,last,k,l 01075 double precision c,f,g,r,s 01076 ! 01077 ! radix=radix floating-point arithmetic de la machine 01078 ! 01079 radix=2.d0 01080 sqrdx=radix**2 01081 l=1 01082 k=n 01083 01084 100 do 2 j=k,1,-1 01085 r=0.d0 01086 do 1 i=1,k 01087 if (i.ne.j) then 01088 r=r+dabs(a(j,i)) 01089 endif 01090 1 continue 01091 if (r.eq.0.d0) then 01092 d(k)=j 01093 if (j.ne.k) then 01094 do 200 i=1,k 01095 f=a(i,j) 01096 a(i,j)=a(i,k) 01097 a(i,k)=f 01098 200 continue 01099 do 201 i=l,k 01100 f=a(j,i) 01101 a(j,i)=a(k,i) 01102 a(k,i)=f 01103 201 continue 01104 endif 01105 k=k-1 01106 goto 2 01107 endif 01108 2 continue 01109 01110 101 do 4 j=l,k 01111 c=0.d0 01112 do 3 i=l,k 01113 if (i.ne.j) then 01114 c=c+dabs(a(i,j)) 01115 endif 01116 3 continue 01117 if (c.eq.0.d0) then 01118 d(l)=j 01119 if (j.ne.l) then 01120 do 202 i=1,k 01121 f=a(i,j) 01122 a(i,j)=a(i,l) 01123 a(i,l)=f 01124 202 continue 01125 do 203 i=l,k 01126 f=a(j,i) 01127 a(j,i)=a(l,i) 01128 a(l,i)=f 01129 203 continue 01130 endif 01131 l=l+1 01132 goto 4 01133 endif 01134 4 continue 01135 01136 linf=l 01137 lsup=k 01138 do 5 i=l,k 01139 d(i)=1.d0 01140 5 continue 01141 01142 102 continue 01143 last=1 01144 do 9 i=l,k 01145 c=0.d0 01146 r=0.d0 01147 ! 01148 ! calcul des normes des colonnes et des lignes 01149 ! 01150 do 6 j=l,k 01151 if (j.ne.i) then 01152 c=c+dabs(a(j,i)) 01153 r=r+dabs(a(i,j)) 01154 endif 01155 6 continue 01156 g=r/radix 01157 f=1.d0 01158 s=c+r 01159 204 if (c.lt.g) then 01160 f=f*radix 01161 c=c*sqrdx 01162 goto 204 01163 endif 01164 g=r*radix 01165 205 if (c.ge.g) then 01166 f=f/radix 01167 c=c/sqrdx 01168 goto 205 01169 endif 01170 if (((c+r)/f).lt.(0.95d0*s)) then 01171 last=0 01172 g=1.d0/f 01173 d(i)=d(i)*f 01174 do 7 j=l,n 01175 a(i,j)=a(i,j)*g 01176 7 continue 01177 do 8 j=1,k 01178 a(j,i)=a(j,i)*f 01179 8 continue 01180 endif 01181 9 continue 01182 if (last.eq.0) goto 102 01183 return 01184 end subroutine balanc 01185 01186 01187 !******************************************************************** 01188 01189 subroutine elmhes(nm,n,linf,lsup,a,intc) 01190 ! 01191 ! reduction de la matrice par la methode d'elimination 01192 ! sous la forme Hessenberg (triangulaire superieure) 01193 ! 01194 integer n,nm,linf,lsup,intc(nm) 01195 double precision a(nm,nm) 01196 ! 01197 ! matrice de dimension n x n 01198 ! 01199 integer i,j,m 01200 double precision x,y 01201 01202 do 17 m=linf+1,lsup-1 01203 x=0.d0 01204 i=m 01205 ! 01206 ! recherche du pivot 01207 ! 01208 do 11 j=m,lsup 01209 if (dabs(a(j,m-1)).gt.dabs(x)) then 01210 x=a(j,m-1) 01211 i=j 01212 endif 01213 11 continue 01214 ! 01215 intc(m)=i 01216 ! 01217 ! echange des colonnes et des lignes 01218 ! 01219 if (i.ne.m) then 01220 do 12 j=m-1,n 01221 y=a(i,j) 01222 a(i,j)=a(m,j) 01223 a(m,j)=y 01224 12 continue 01225 do 13 j=1,lsup 01226 y=a(j,i) 01227 a(j,i)=a(j,m) 01228 a(j,m)=y 01229 13 continue 01230 endif 01231 ! 01232 ! elimination 01233 ! 01234 if (x.ne.0.d0) then 01235 do 16 i=m+1,lsup 01236 y=a(i,m-1) 01237 if (y.ne.0.d0) then 01238 y=y/x 01239 a(i,m-1)=y 01240 do 14 j=m,n 01241 a(i,j)=a(i,j)-y*a(m,j) 01242 14 continue 01243 do 15 j=1,lsup 01244 a(j,m)=a(j,m)+y*a(j,i) 01245 15 continue 01246 endif 01247 16 continue 01248 endif 01249 17 continue 01250 01251 return 01252 end subroutine elmhes 01253 01254 01255 !******************************************************************** 01256 01257 subroutine eltran(nm,n,linf,lsup,a,intc,vecs) 01258 01259 integer i,j,k 01260 integer n,nm,linf,lsup,intc(nm) 01261 double precision a(nm,nm),vecs(nm,nm) 01262 01263 do 2 i=1,n 01264 do 1 j=1,n 01265 vecs(i,j)=0.d0 01266 1 continue 01267 vecs(i,i)=1.d0 01268 2 continue 01269 01270 do 5 i=lsup-1,linf+1,-1 01271 j=intc(i) 01272 do 3 k=i+1,lsup 01273 vecs(k,i)=a(k,i-1) 01274 3 continue 01275 if (i.ne.j) then 01276 do 4 k=i,lsup 01277 vecs(i,k)=vecs(j,k) 01278 vecs(j,k)=0.d0 01279 4 continue 01280 vecs(j,i)=1.d0 01281 endif 01282 5 continue 01283 01284 return 01285 end subroutine eltran 01286 01287 !******************************************************************** 01288 01289 subroutine hqr2(nm,n,linf,lsup,a,vecs,wr,wi) 01290 01291 integer n,nm,linf,lsup 01292 double precision a(nm,nm),wr(nm),wi(nm),vecs(nm,nm) 01293 01294 01295 integer i,its,j,k,l,m,na,en 01296 double precision anorm,p,q,r,s,t,u,v,w,x,y,z,ra 01297 double precision sa,vr,vi,eps 01298 01299 do 1 i=1,linf-1 01300 wr(i)=a(i,i) 01301 wi(i)=0.d0 01302 1 continue 01303 do 2 i=lsup+1,n 01304 wr(i)=a(i,i) 01305 wi(i)=0.d0 01306 2 continue 01307 01308 eps=2.220D-16 01309 en=lsup 01310 t=0.d0 01311 100 if (en.lt.linf) goto 200 01312 its=0 01313 na=en-1 01314 101 do 99 l=en,linf+1,-1 01315 s=dabs(a(l-1,l-1))+dabs(a(l,l)) 01316 if ((dabs(a(l,l-1))+s).eq.s) goto 102 01317 99 continue 01318 l=linf 01319 102 x=a(en,en) 01320 if (l.eq.en) goto 201 01321 y=a(na,na) 01322 w=a(en,na)*a(na,en) 01323 if (l.eq.na) goto 202 01324 if (its.eq.30) then 01325 write(0,*) 'nombre max iter atteint' 01326 return 01327 endif 01328 if (its.eq.10.or.its.eq.20) then 01329 t=t+x 01330 do 103 i=linf,en 01331 a(i,i)=a(i,i)-x 01332 103 continue 01333 s=dabs(a(en,na))+dabs(a(na,en-2)) 01334 x=0.75d0*s 01335 y=x 01336 w=-0.4375d0*s*s 01337 endif 01338 its=its+1 01339 do 104 m=en-2,l,-1 01340 z=a(m,m) 01341 r=x-z 01342 s=y-z 01343 p=(r*s-w)/a(m+1,m)+a(m,m+1) 01344 q=a(m+1,m+1)-z-r-s 01345 r=a(m+2,m+1) 01346 s=dabs(p)+dabs(q)+dabs(r) 01347 p=p/s 01348 q=q/s 01349 r=r/s 01350 if (m.eq.l) goto 105 01351 u=dabs(a(m,m-1))*(dabs(q)+dabs(r)) 01352 v=dabs(p)*(dabs(a(m-1,m-1))+dabs(z)+dabs(a(m+1,m+1))) 01353 if ((u+v).eq.v) goto 105 01354 104 continue 01355 105 do 106 i=m+2,en 01356 a(i,i-2)=0.d0 01357 if (i.ne.m+2) a(i,i-3)=0.d0 01358 106 continue 01359 do 110 k=m,na 01360 if (k.ne.m) then 01361 p=a(k,k-1) 01362 q=a(k+1,k-1) 01363 r=0.d0 01364 if (k.ne.na) r=a(k+2,k-1) 01365 x=dabs(p)+dabs(q)+dabs(r) 01366 if (x.eq.0) goto 110 01367 p=p/x 01368 q=q/x 01369 r=r/x 01370 endif 01371 s=DSQRT(p**2+q**2+r**2) 01372 if (p.lt.0.d0) s=-s 01373 if (k.ne.m) then 01374 a(k,k-1)=-s*x 01375 else 01376 if (l.ne.m) a(k,k-1)=-a(k,k-1) 01377 endif 01378 p=p+s 01379 x=p/s 01380 y=q/s 01381 z=r/s 01382 q=q/p 01383 r=r/p 01384 do 107 j=k,n 01385 p=a(k,j)+q*a(k+1,j) 01386 if (k.ne.na) then 01387 p=p+r*a(k+2,j) 01388 a(k+2,j)=a(k+2,j)-p*z 01389 endif 01390 a(k+1,j)=a(k+1,j)-p*y 01391 a(k,j)=a(k,j)-p*x 01392 107 continue 01393 do 108 i=1,min(en,k+3) 01394 p=x*a(i,k)+y*a(i,k+1) 01395 if (k.ne.na) then 01396 p=p+z*a(i,k+2) 01397 a(i,k+2)=a(i,k+2)-p*r 01398 endif 01399 a(i,k+1)=a(i,k+1)-p*q 01400 a(i,k)=a(i,k)-p 01401 108 continue 01402 do 109 i=linf,lsup 01403 p=x*vecs(i,k)+y*vecs(i,k+1) 01404 if (k.ne.na) then 01405 p=p+z*vecs(i,k+2) 01406 vecs(i,k+2)=vecs(i,k+2)-p*r 01407 endif 01408 vecs(i,k+1)=vecs(i,k+1)-p*q 01409 vecs(i,k)=vecs(i,k)-p 01410 109 continue 01411 110 continue 01412 goto 101 01413 01414 201 wr(en)=x+t 01415 a(en,en)=wr(en) 01416 wi(en)=0.d0 01417 en=na 01418 goto 100 01419 01420 202 p=(y-x)/2.d0 01421 q=p**2+w 01422 z=DSQRT(dabs(q)) 01423 x=x+t 01424 a(en,en)=x 01425 a(na,na)=y+t 01426 if (q.gt.0.d0) then 01427 if (p.lt.0.d0) then 01428 z=p-z 01429 else 01430 z=p+z 01431 endif 01432 wr(na)=x+z 01433 wr(en)=x-w/z 01434 s=wr(en) 01435 wi(na)=0.d0 01436 wi(en)=wi(na) 01437 x=a(en,na) 01438 r=DSQRT(x**2+z**2) 01439 p=x/r 01440 q=z/r 01441 do 111 j=na,n 01442 z=a(na,j) 01443 a(na,j)=q*z+p*a(en,j) 01444 a(en,j)=q*a(en,j)-p*z 01445 111 continue 01446 do 112 i=1,en 01447 z=a(i,na) 01448 a(i,na)=q*z+p*a(i,en) 01449 a(i,en)=q*a(i,en)-p*z 01450 112 continue 01451 do 113 i=linf,lsup 01452 z=vecs(i,na) 01453 vecs(i,na)=q*z+p*vecs(i,en) 01454 vecs(i,en)=q*vecs(i,en)-p*z 01455 113 continue 01456 else 01457 wr(na)=x+p 01458 wr(en)=wr(na) 01459 wi(na)=z 01460 wi(en)=-z 01461 endif 01462 en=en-2 01463 goto 100 01464 01465 200 anorm=0.d0 01466 k=1 01467 do 115 i=1,n 01468 do 114 j=k,n 01469 anorm=anorm+dabs(a(i,j)) 01470 114 continue 01471 k=i 01472 115 continue 01473 01474 do 120 en=n,1,-1 01475 p=wr(en) 01476 q=wi(en) 01477 na=en-1 01478 if (q.eq.0.d0) then 01479 m=en 01480 a(en,en)=1.d0 01481 do 117 i=na,1,-1 01482 w=a(i,i)-p 01483 r=a(i,en) 01484 do 116 j=m,na 01485 r=r+a(i,j)*a(j,en) 01486 116 continue 01487 if (wi(i).lt.0.d0) then 01488 z=w 01489 s=r 01490 else 01491 m=i 01492 if (wi(i).eq.0.d0) then 01493 if (w.ne.0) then 01494 a(i,en)=-r/w 01495 else 01496 a(i,en)=-r/eps*anorm 01497 endif 01498 else 01499 x=a(i,i+1) 01500 y=a(i+1,i) 01501 q=(wr(i)-p)**2+wi(i)**2 01502 a(i,en)=(x*s-z*r)/q 01503 t=a(i,en) 01504 if (dabs(x).gt.dabs(z)) then 01505 a(i+1,en)=(-r-w*t)/x 01506 else 01507 a(i+1,en)=(-s-y*t)/z 01508 endif 01509 endif 01510 endif 01511 117 continue 01512 else 01513 if (q.lt.0.d0) then 01514 m=na 01515 if (dabs(a(en,na)).gt.dabs(a(na,en))) then 01516 a(na,na)=-(a(en,en)-p)/a(en,na) 01517 a(na,en)=-q/a(en,na) 01518 else 01519 call cdiv(-a(na,en),0.d0,a(na,na)-p,q,a(na,na),a(na,en)) 01520 endif 01521 a(en,na)=1.d0 01522 a(en,en)=0.d0 01523 do 119 i=na-1,1,-1 01524 w=a(i,i)-p 01525 ra=a(i,en) 01526 sa=0.d0 01527 do 118 j=m,na 01528 ra=ra+a(i,j)*a(j,na) 01529 sa=sa+a(i,j)*a(j,en) 01530 118 continue 01531 if (wi(i).lt.0.d0) then 01532 z=w 01533 r=ra 01534 s=sa 01535 else 01536 m=i 01537 if (wi(i).eq.0.d0) then 01538 call cdiv(-ra,-sa,w,q,a(i,na),a(i,en)) 01539 else 01540 x=a(i,i+1) 01541 y=a(i+1,i) 01542 vr=(wr(i)-p)**2+wi(i)**2-q**2 01543 vi=(wr(i)-p)*2.d0*q 01544 if ((vr.eq.0.d0).and.(vi.eq.0.d0)) then 01545 vr=eps*anorm*(dabs(w)+dabs(q)+dabs(x)+dabs(y)+dabs(z)) 01546 call cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi,a(i,na),a(i,en)) 01547 endif 01548 if (dabs(x).gt.(dabs(z)+dabs(q))) then 01549 a(i+1,na)=(-ra-w*a(i,na)+q*a(i,en))/x 01550 a(i+1,en)=(-sa-w*a(i,en)-q*a(i,na))/x 01551 else 01552 call cdiv(-r-y*a(i,na),-s-y*a(i,en),z,q,a(i+1,na),a(i+1,en)) 01553 endif 01554 endif 01555 endif 01556 119 continue 01557 endif 01558 endif 01559 120 continue 01560 01561 do 122 i=1,linf-1 01562 do 121 j=i+1,n 01563 vecs(i,j)=a(i,j) 01564 121 continue 01565 122 continue 01566 do 124 i=lsup+1,n 01567 do 123 j=i+1,n 01568 vecs(i,j)=a(i,j) 01569 123 continue 01570 124 continue 01571 01572 do 129 j=n,linf,-1 01573 if (j.le.lsup) then 01574 m=j 01575 else 01576 m=lsup 01577 endif 01578 l=j-1 01579 if (wi(j).lt.0.d0) then 01580 do 126 i=linf,lsup 01581 y=0.d0 01582 z=0.d0 01583 do 125 k=linf,m 01584 y=y+vecs(i,k)*vecs(k,l) 01585 z=z+vecs(i,k)*vecs(k,j) 01586 125 continue 01587 vecs(i,l)=y 01588 vecs(i,j)=z 01589 126 continue 01590 else 01591 if (wi(j).eq.0.d0) then 01592 do 128 i=linf,lsup 01593 z=0.d0 01594 do 127 k=linf,m 01595 z=z+vecs(i,k)*a(k,j) 01596 127 continue 01597 vecs(i,j)=z 01598 128 continue 01599 endif 01600 endif 01601 129 continue 01602 01603 return 01604 end subroutine hqr2 01605 01606 !******************************************************************** 01607 01608 subroutine balbak(nm,n,linf,lsup,d,z) 01609 01610 integer n,nm,linf,lsup 01611 double precision d(nm),z(nm,nm) 01612 01613 integer i,j,k 01614 double precision s 01615 01616 do 2 i=linf,lsup 01617 s=d(i) 01618 do 1 j=1,n 01619 z(i,j)=s*z(i,j) 01620 1 continue 01621 2 continue 01622 01623 do 4 i=linf-1,1,-1 01624 k=d(i) 01625 if (k.ne.i) then 01626 do 3 j=1,n 01627 s=z(i,j) 01628 z(i,j)=z(k,j) 01629 z(k,j)=s 01630 3 continue 01631 endif 01632 4 continue 01633 do 6 i=lsup+1,n 01634 k=d(i) 01635 if (k.ne.i) then 01636 do 5 j=1,n 01637 s=z(i,j) 01638 z(i,j)=z(k,j) 01639 z(k,j)=s 01640 5 continue 01641 endif 01642 6 continue 01643 01644 return 01645 end subroutine balbak 01646 01647 !******************************************************************** 01648 01649 subroutine cdiv(xr,xi,yr,yi,zr,zi) 01650 01651 double precision xr,xi,yr,yi,zr,zi,h 01652 01653 if (yr.eq.0.d0.or.yi.eq.0.d0) then 01654 write(0,*) 'pb dans cdiv' 01655 return 01656 endif 01657 01658 if (dabs(yr).gt.dabs(yi)) then 01659 h=yi/yr 01660 yr=h*yi+yr 01661 zr=(xr+h*xi)/yr 01662 zi=(xi-h*xr)/yr 01663 else 01664 h=yr/yi 01665 yi=h*yr+yi 01666 zr=(h*xr+xi)/yi 01667 zi=(h*xi-xr)/yi 01668 endif 01669 return 01670 end subroutine cdiv 01671 01672 01673 !******************************************************************** 01674 01675 end module m_vpropre 01676 01677 !------------------------------------------------------------------------------ 01678 ! 01679 ! VERSION : 1.0 01680 ! 01681 ! MODULE: tools.f90 function RNRNOR 01702 DOUBLE PRECISION FUNCTION RNRNOR() 01703 !C 01704 !C RNRNOR generates normal random numbers with zero mean and unit 01705 !C standard deviation, often denoted N(0,1) 01706 !C 01707 INTEGER J, N, TN 01708 DOUBLE PRECISION TWOPIS, AA, B, C, XDN 01709 PARAMETER ( N = 64, TN = 2*N, TWOPIS = TN/2.506628274631000D0 ) 01710 PARAMETER ( XDN = 0.3601015713011893D0, B = 0.4878991777603940D0 ) 01711 PARAMETER ( AA = 12.37586029917064D0, C = 12.67705807886560D0 ) 01712 DOUBLE PRECISION XT, XX, Y, RNRUNI 01713 DOUBLE PRECISION X(0:N) 01714 SAVE X 01715 DATA ( X(J), J = 0, 31 ) / & 01716 0.3409450287039653D+00, 0.4573145918669259D+00, & 01717 0.5397792816116612D+00, 0.6062426796530441D+00, & 01718 0.6631690627645207D+00, 0.7136974590560222D+00, & 01719 0.7596124749339174D+00, 0.8020356003555283D+00, & 01720 0.8417226679789527D+00, 0.8792102232083114D+00, & 01721 0.9148948043867484D+00, 0.9490791137530882D+00, & 01722 0.9820004812398864D+00, 0.1013849238029940D+01, & 01723 0.1044781036740172D+01, 0.1074925382028552D+01, & 01724 0.1104391702268125D+01, 0.1133273776243940D+01, & 01725 0.1161653030133931D+01, 0.1189601040838737D+01, & 01726 0.1217181470700870D+01, 0.1244451587898246D+01, & 01727 0.1271463480572119D+01, 0.1298265041883197D+01, & 01728 0.1324900782180860D+01, 0.1351412509933371D+01, & 01729 0.1377839912870011D+01, 0.1404221063559975D+01, & 01730 0.1430592868502691D+01, 0.1456991476137671D+01, & 01731 0.1483452656603219D+01, 0.1510012164318519D+01 / 01732 DATA ( X(J), J = 32, 64 ) / & 01733 0.1536706093359520D+01, 0.1563571235037691D+01, & 01734 0.1590645447014253D+01, 0.1617968043674446D+01, & 01735 0.1645580218369081D+01, 0.1673525509567038D+01, & 01736 0.1701850325062740D+01, 0.1730604541317782D+01, & 01737 0.1759842199038300D+01, 0.1789622321566574D+01, & 01738 0.1820009890130691D+01, 0.1851077020230275D+01, & 01739 0.1882904397592872D+01, 0.1915583051943031D+01, & 01740 0.1949216574916360D+01, 0.1983923928905685D+01, & 01741 0.2019843052906235D+01, 0.2057135559990095D+01, & 01742 0.2095992956249391D+01, 0.2136645022544389D+01, & 01743 0.2179371340398135D+01, 0.2224517507216017D+01, & 01744 0.2272518554850147D+01, 0.2323933820094302D+01, & 01745 0.2379500774082828D+01, 0.2440221797979943D+01, & 01746 0.2507511701865317D+01, 0.2583465835225429D+01, & 01747 0.2671391590320836D+01,4*0.2776994269662875D+01 / 01748 Y = RNRUNI() 01749 J = MOD( INT( TN*RNRUNI() ), N ) 01750 XT = X(J+1) 01751 RNRNOR = ( Y + Y - 1 )*XT 01752 IF ( DABS(RNRNOR) .GT. X(J) ) THEN 01753 XX = B*( XT - DABS(RNRNOR) )/( XT - X(J) ) 01754 Y = RNRUNI() 01755 IF ( Y .GT. C - AA*DEXP( -XX**2/2 ) ) THEN 01756 RNRNOR = SIGN( XX, RNRNOR ) 01757 ELSE 01758 IF ( DEXP(-XT**2/2)+Y/(TWOPIS*XT).GT.DEXP(-RNRNOR**2/2) ) THEN 01759 10 XX = XDN*DLOG( RNRUNI() ) 01760 IF ( -2*DLOG( RNRUNI() ) .LE. XX**2 ) GO TO 10 01761 RNRNOR = SIGN( X(N) - XX, RNRNOR ) 01762 END IF 01763 END IF 01764 END IF 01765 END FUNCTION RNRNOR 01766 01767 !------------------------------------------------------------------------------ 01768 ! 01769 ! VERSION : 1.0 01770 ! 01771 ! MODULE: tools.f90 function RNRUNI 01791 DOUBLE PRECISION FUNCTION RNRUNI() 01792 !C 01793 !C Uniform (0, 1) random number generator 01794 !C 01795 INTEGER A12, A13, A21, A23, P12, P13, P21, P23 01796 INTEGER Q12, Q13, Q21, Q23, R12, R13, R21, R23 01797 INTEGER X10, X11, X12, X20, X21, X22, Z, M1, M2, H 01798 DOUBLE PRECISION INVMP1 01799 PARAMETER ( M1 = 2147483647, M2 = 2145483479 ) 01800 PARAMETER ( A12 = 63308, Q12 = 33921, R12 = 12979 ) 01801 PARAMETER ( A13 = -183326, Q13 = 11714, R13 = 2883 ) 01802 PARAMETER ( A21 = 86098, Q21 = 24919, R21 = 7417 ) 01803 PARAMETER ( A23 = -539608, Q23 = 3976, R23 = 2071 ) 01804 PARAMETER ( INVMP1 = 4.656612873077392578125D-10 ) 01805 01806 SAVE X10, X11, X12, X20, X21, X22 01807 DATA X10, X11, X12, X20, X21, X22 & 01808 / 15485857, 17329489, 36312197, 55911127, 75906931, 96210113 / 01809 !C 01810 !C Component 1 01811 !C 01812 H = X10/Q13 01813 P13 = -A13*( X10 - H*Q13 ) - H*R13 01814 H = X11/Q12 01815 P12 = A12*( X11 - H*Q12 ) - H*R12 01816 IF ( P13 .LT. 0 ) P13 = P13 + M1 01817 IF ( P12 .LT. 0 ) P12 = P12 + M1 01818 X10 = X11 01819 X11 = X12 01820 X12 = P12 - P13 01821 IF ( X12 .LT. 0 ) X12 = X12 + M1 01822 !C 01823 !C Component 2 01824 !C 01825 H = X20/Q23 01826 P23 = -A23*( X20 - H*Q23 ) - H*R23 01827 H = X22/Q21 01828 P21 = A21*( X22 - H*Q21 ) - H*R21 01829 IF ( P23 .LT. 0 ) P23 = P23 + M2 01830 IF ( P21 .LT. 0 ) P21 = P21 + M2 01831 X20 = X21 01832 X21 = X22 01833 X22 = P21 - P23 01834 IF ( X22 .LT. 0 ) X22 = X22 + M2 01835 !C 01836 !C Combination 01837 !C 01838 Z = X12 - X22 01839 IF ( Z .LE. 0 ) Z = Z + M1 01840 RNRUNI = Z*INVMP1 01841 END FUNCTION RNRUNI 01842 01843 01844
1.7.3