N.I.M.R.O.D.

tools.f90

Go to the documentation of this file.
```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
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
01073
01074            integer i,j,last,k,l
01075            double precision c,f,g,r,s
01076            !
01078            !
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
01157            f=1.d0
01158            s=c+r
01159 204        if (c.lt.g) then
01161                c=c*sqrdx
01162                goto 204
01163            endif
01165 205        if (c.ge.g) then
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
```