DOUG 0.2

pcg.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !!----------------------------------
00023 !! Preconditioned Conjugate Gradient
00024 !!----------------------------------
00025 module pcg_mod
00026 
00027   use ConvInf_mod
00028   use SpMtx_mods
00029   use Mesh_class
00030   use globals
00031   use subsolvers
00032   use Preconditioner_mod
00033   use Distribution_mod
00034 
00035   implicit none
00036 
00037 #include<doug_config.h>
00038 
00039 #ifdef D_COMPLEX
00040 #define float complex
00041 #else
00042 #define float real
00043 #endif
00044 
00045   private :: &
00046        preconditioner, &
00047        msolve
00048 
00049   ! profiling
00050   real(kind=rk), private, save :: time_preconditioner = 0
00051 
00052 contains
00053 
00054   subroutine prec2Level(prepare,A,CP,sol,rhs,res)
00055     use CoarseAllgathers
00056 
00057     logical, intent(in) :: prepare
00058     type(SpMtx)  :: A   !< sparse system matrix
00059     type(CoarsePreconditioner),intent(inout) :: CP
00060     real(kind=rk),dimension(:),pointer :: sol !< solution
00061     real(kind=rk),dimension(:),pointer :: rhs !< right hand side
00062     real(kind=rk),dimension(:),pointer :: res !< residual vector
00063 
00064     real(kind=rk),dimension(:),pointer,save :: csol,crhs,clrhs,tmpsol2,tmpsol
00065     integer :: ol
00066 
00067     ol=max(sctls%overlap,sctls%smoothers)
00068 
00069     if (prepare) then
00070       if (sctls%verbose>4) write(stream,*) "Preparing 2. level"
00071       call prec2Level_prepare()
00072     else
00073       if (sctls%verbose>4) write(stream,*) "Solving 2. level"
00074       call prec2Level_solve()
00075     end if
00076   
00077   contains
00078     subroutine prec2Level_exchangeMatrix()
00079       integer :: ol
00080       logical :: add
00081 
00082       if (sctls%verbose>4) write(stream,*) "Exchanging coarse matrix"
00083 
00084       ol=max(sctls%overlap,sctls%smoothers)
00085       if (ol==0) then
00086         add=.false.
00087       else
00088         add=.true.
00089       endif
00090 
00091       call AllSendCoarseMtx(CP%cdat%LAC,CP%AC,CP%cdat%lg_cfmap,&
00092            CP%cdat%ngfc,CP%cdat%nprocs,CP%cdat%send)
00093       call AllRecvCoarseMtx(CP%AC,CP%cdat%send,add=add) ! Recieve it
00094 
00095     end subroutine prec2Level_exchangeMatrix
00096 
00097     subroutine prec2Level_prepare()
00098       if (.not.CP%ready) then
00099         if (CP%cdat%active) then
00100           ! First iteration - send matrix
00101           call prec2Level_exchangeMatrix()
00102         end if
00103 
00104         ! after coarse matrix is received allocate coarse vectors
00105         if (associated(crhs)) then
00106           if (size(crhs)/=CP%AC%ncols) then
00107             deallocate(csol)
00108             deallocate(crhs)
00109             allocate(crhs(CP%AC%ncols))
00110             allocate(csol(CP%AC%nrows))
00111           endif
00112         else
00113           allocate(crhs(CP%AC%ncols))
00114           allocate(csol(CP%AC%nrows))
00115         endif
00116         ! allocate memory for vector
00117         if (CP%cdat%active) then
00118           allocate(clrhs(CP%cdat%nlfc))
00119         else
00120           allocate(clrhs(CP%cdat_vec%nlfc))
00121         end if
00122         CP%ready = .true.
00123       end if
00124 
00125       if (.not.associated(tmpsol)) then
00126         !allocate(tmpsol(A%nrows))
00127         allocate(tmpsol(size(rhs)))
00128       endif
00129 
00130       if (CP%cdat%active) then
00131         if (sctls%verbose>6) write(stream,*) "Restricting into local coarse vector", size(clrhs)
00132         ! Send coarse vector
00133         call SpMtx_Ax(clrhs,CP%R,rhs,dozero=.true.) ! restrict <RA>
00134         if (CP%cdat_vec%active) then
00135           call AllSendCoarseVector(clrhs,CP%cdat_vec%nprocs,CP%cdat_vec%cdisps,&
00136                CP%cdat_vec%send)
00137         else
00138           call AllSendCoarseVector(clrhs,CP%cdat%nprocs,CP%cdat%cdisps,&
00139                CP%cdat%send)
00140         endif
00141       end if ! CP%cdat%active
00142     end subroutine prec2Level_prepare
00143 
00144     subroutine prec2Level_solve()
00145       if (.not.CP%cdat%active) then ! 1 processor case
00146         if (sctls%method>1.and.sctls%method/=5) then ! multiplicative Schwarz
00147           call SpMtx_Ax(crhs,CP%R,res,dozero=.true.) ! restriction
00148         else
00149           call SpMtx_Ax(crhs,CP%R,rhs,dozero=.true.) ! restriction
00150         endif
00151       end if
00152 
00153       !csol=0.0_rk
00154       if (CP%cdat%active) then
00155         ! Recieve the vector for solve
00156         if (CP%cdat_vec%active) then
00157           call AllRecvCoarseVector(crhs,CP%cdat_vec%nprocs,&
00158                CP%cdat_vec%cdisps,CP%cdat_vec%glg_cfmap,CP%cdat_vec%send)
00159         else
00160           call AllRecvCoarseVector(crhs,CP%cdat%nprocs,&
00161                CP%cdat%cdisps,CP%cdat%glg_cfmap,CP%cdat%send)
00162         endif
00163         !call MPI_BARRIER(MPI_COMM_WORLD,i)
00164       end if
00165 
00166       if (CP%AC%subsolve_id<=0) then
00167         write (stream,*) &
00168              'factorising coarse matrix of size',CP%AC%nrows, &
00169              ' and nnz:',CP%AC%nnz
00170       end if
00171 
00172       ! Coarse solve
00173       call sparse_singlesolve(CP%AC%subsolve_id,csol,crhs,&
00174            nfreds=CP%AC%nrows, &
00175            nnz=CP%AC%nnz,        &
00176            indi=CP%AC%indi,      &
00177            indj=CP%AC%indj,      &
00178            val=CP%AC%val)
00179 
00180       if (CP%cdat_vec%active) then
00181         call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.)
00182         call SpMtx_Ax(tmpsol,CP%R,clrhs,dozero=.true.,transp=.true.)
00183 
00184       elseif (CP%cdat%active) then
00185         call Vect_remap(csol,clrhs,CP%cdat%gl_cfmap,dozero=.true.)
00186         call SpMtx_Ax(tmpsol,CP%R,clrhs,dozero=.true.,transp=.true.)
00187       else
00188         call SpMtx_Ax(tmpsol,CP%R,csol,dozero=.true.,transp=.true.) ! interpolation
00189       endif
00190 
00191       if (sctls%method==1) then
00192         !call Print_Glob_Vect(sol,M,'sol===',chk_endind=M%ninner)
00193         sol=sol+tmpsol
00194       elseif (sctls%method==3) then ! fully multiplicative Schwarz
00195         sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
00196         ! calculate the residual:
00197         call SpMtx_Ax(res,A,sol,dozero=.true.) ! 
00198         res=rhs-res
00199       elseif (sctls%method==2.and.ol>0) then ! fully multiplicative Schwarz
00200         sol(1:A%nrows)=tmpsol(1:A%nrows)
00201         ! calculate the residual:
00202         call SpMtx_Ax(res,A,sol,dozero=.true.) ! 
00203         res=rhs-res
00204       endif
00205       if (((ol==0.and.sctls%method==2).or.sctls%method==5).and.sctls%levels>1) then 
00206         ! multiplicative on fine level, additive with coarse level: 
00207         sol(1:A%nrows)=sol(1:A%nrows)+tmpsol(1:A%nrows)
00208       endif
00209     end subroutine prec2Level_solve
00210 
00211   end subroutine prec2Level
00212 
00213   !-------------------------------
00216   subroutine preconditioner(sol,A,rhs,M,finePrec,coarsePrec,&
00217                A_ghost,bugtrack_)
00218     use CoarseAllgathers
00219     use CoarseMtx_mod
00220     use Vect_mod
00221     implicit none
00222     real(kind=rk),dimension(:),pointer :: sol !< solution
00223     type(SpMtx)                        :: A   !< sparse system matrix
00224     real(kind=rk),dimension(:),pointer :: rhs !< right hand side
00225     type(Mesh),intent(in)              :: M   !< Mesh
00226     type(FinePreconditioner),intent(inout) :: finePrec !< fine level preconditioner
00227     type(CoarsePreconditioner),intent(inout) :: coarsePrec !< coarse level preconditioner
00228     real(kind=rk),dimension(:),pointer :: res !< residual vector, allocated
00229                                               !! here for multiplicative Schwarz
00230     type(SpMtx),optional               :: A_ghost  !< matr@interf.
00231     logical,optional                   :: bugtrack_
00232     ! ----- local: ------
00233     integer :: i
00234     logical :: add,bugtrack
00235     real(kind=rk) :: t1
00236 
00237     if (sctls%verbose>4) write(stream,*) "Applying preconditioner"
00238     t1 = MPI_WTime()
00239 
00240     if (present(bugtrack_)) then
00241       bugtrack=bugtrack_
00242     else
00243       bugtrack=.false.
00244     endif
00245     ! ----------------------------
00246 
00247     if (sctls%method==0) then
00248       sol=rhs
00249       return
00250     endif
00251 
00252     sol=0.0_rk
00253     if (sctls%method>1) then ! For multiplicative Schwarz method...:
00254       allocate(res(size(rhs)))
00255     endif
00256       
00257     if (sctls%levels>1) then
00258       call prec2Level(.true.,A,coarsePrec,sol,rhs,res)
00259     end if
00260 
00261     ! first level prec
00262     call FinePreconditioner_apply(finePrec,sol,rhs)
00263 
00264     if (sctls%levels>1) then
00265       call prec2Level(.false.,A,coarsePrec,sol,rhs,res)
00266     end if
00267 
00268     time_preconditioner = time_preconditioner + MPI_WTime()-t1
00269 
00270   end subroutine preconditioner
00271 
00272   !--------------------------
00273   ! Solve
00274   !--------------------------
00275   subroutine msolve (m, r, z)
00276     implicit none
00277     real(kind=rk), dimension(:),     intent(in) :: m, r
00278     real(kind=rk), dimension(:), intent(in out) :: z
00279     integer :: i
00280 
00281     do i = 1, size(r)
00282        z(i) = r(i) !/ m(i)
00283     end do
00284   end subroutine msolve
00285 
00286   !--------------------------
00289   subroutine pcg_weigs (D,x,finePrec,coarsePrec,it,cond_num,tol_,maxit_, &
00290        x0_,solinf,resvects_)
00291     use CoarseAllgathers
00292 
00293     implicit none
00294     
00295     type(Distribution),intent(inout) :: D
00296     float(kind=rk),dimension(:),pointer :: x !< Solution
00297     type(FinePreconditioner),intent(inout)   :: finePrec !< fine level preconditioner
00298     type(CoarsePreconditioner),intent(inout) :: coarsePrec !< coarse level preconditioner
00299 
00300     integer,intent(out) :: it
00301     real(kind=rk),intent(out) :: cond_num
00302     
00303     ! Optional arguments
00304     real(kind=rk), intent(in), optional                :: tol_ !< Tolerance of the method
00305     integer,       intent(in), optional                :: maxit_ !< Max number of iterations
00306     float(kind=rk), dimension(:), intent(in), optional :: x0_ !< Initial guess
00307     type(ConvInf),            intent(in out), optional :: solinf !< Solution statistics
00308     logical,                      intent(in), optional :: resvects_ !< Fill in the 'resvect' or not
00309     
00310     real(kind=rk) :: tol   ! Tolerance
00311     integer       :: maxit ! Max number of iterations
00312     logical       :: resvects=.false.
00313 
00314     real(kind=rk),dimension(:),pointer,save :: r, p, q, m, z, b_k
00315     real(kind=rk),dimension(:),pointer,save :: ztmp !TODO remove me
00316     real(kind=rk) :: rho_curr, rho_prev, init_norm, res_norm
00317     real(kind=rk) :: tmp,ratio_norm
00318     integer :: nfreds
00319     integer :: ierr
00320     real(kind=rk),dimension(:),pointer,save :: dd,ee,alpha,beta
00321     !logical,parameter :: bugtrack=.true.
00322     logical,parameter :: bugtrack=.false.
00323     ! profiling
00324     real(8) :: time_iterations, t1
00325 
00326     write(stream,'(/a)') 'Preconditioned conjugate gradient:'
00327 
00328     if (size(D%rhs) /= size(x)) &
00329          call DOUG_abort('[pcg] : SEVERE : size(b) /= size(x)',-1)
00330 
00331     nfreds=size(D%rhs)
00332     if (.not.associated(r)) then
00333       allocate(r(nfreds))
00334       allocate(p(nfreds))
00335       allocate(q(nfreds))
00336       allocate(m(nfreds))
00337       allocate(z(nfreds))
00338       allocate(b_k(nfreds))
00339     endif
00340 
00341     tol = sctls%solve_tolerance
00342     if (present(tol_)) &
00343          tol = tol_
00344 
00345     maxit = sctls%solve_maxiters
00346     if (maxit==-1) maxit=100000
00347     if (present(maxit_)) &
00348          maxit = maxit_
00349 
00350     if (present(resvects_).and.(.not.present(solinf))) &
00351          call DOUG_abort('[pcg] : SEVERE : "resvects_" must be given along'//&
00352          ' with "solinf".',-1)
00353 
00354     ! Allocate convergence info structure
00355     if (present(solinf).and.(.not.present(resvects_))) then
00356        call ConvInf_Init(solinf)
00357     else if (present(solinf).and.&
00358          (present(resvects_).and.(resvects_.eqv.(.true.)))) then
00359        call ConvInf_Init(solinf, maxit)
00360        resvects = .true.
00361     end if
00362 
00363 if (bugtrack)call Print_Glob_Vect(x,D%mesh,'global x===')
00364     call Distribution_pmvm(D,r,x)
00365 
00366     r = D%rhs - r
00367     init_norm = Vect_dot_product(r,r)
00368     if (init_norm == 0.0) &
00369          init_norm = 1.0
00370 
00371 
00372     write(stream,'(a,i6,a,e8.2)') 'maxit = ',maxit,', tol = ',tol
00373     write(stream,'(a,e22.15)') 'init_norm = ', dsqrt(init_norm)
00374 
00375     it = 0
00376     ratio_norm = 1.0_rk
00377     ! arrays for eigenvalue calculation:
00378     if (.not.associated(dd)) then
00379       allocate(dd(maxit))
00380       allocate(ee(maxit))
00381       allocate(alpha(maxit))
00382       allocate(beta(maxit))
00383     endif
00384 
00385     time_iterations = 0.
00386     ! iterations
00387     do while((ratio_norm > tol*tol).and.(it < maxit))
00388       it = it + 1
00389       t1 = MPI_WTime()
00390 
00391 if (bugtrack)call Print_Glob_Vect(r,D%mesh,'global r===',chk_endind=D%mesh%ninner)
00392       call preconditioner(sol=z,          &
00393                             A=D%A,          &
00394                           rhs=r,          &
00395                             M=D%mesh,        &
00396                             finePrec=finePrec, &
00397                             coarsePrec=coarsePrec, &
00398                     A_ghost=D%A_ghost,  &
00399                     bugtrack_=bugtrack)
00400 
00401 if (bugtrack)call Print_Glob_Vect(z,D%mesh,'global bef comm z===',chk_endind=D%mesh%ninner)
00402       if (sctls%method/=0) then
00403         call Distribution_addoverlap(D,z)
00404       endif
00405 
00406 !call Print_Glob_Vect(z,D%mesh,'global aft comm z===',chk_endind=D%mesh%ninner)
00407 !call Print_Glob_Vect(z,D%mesh,'global z===')
00408       ! compute current rho
00409       rho_curr = Vect_dot_product(r,z)
00410 
00411       if (it == 1) then
00412          p = z
00413       else
00414          beta(it) = rho_curr / rho_prev
00415          p = z + beta(it) * p
00416       end if
00417       call Distribution_pmvm(D,q,p)
00418 if (bugtrack)call Print_Glob_Vect(q,D%mesh,'global q===')
00419       ! compute alpha
00420       alpha(it) = rho_curr / Vect_dot_product(p,q)
00421       x = x + alpha(it) * p
00422 if (bugtrack)call Print_Glob_Vect(x,D%mesh,'global x===')
00423       r = r - alpha(it) * q
00424 if (bugtrack)call Print_Glob_Vect(r,D%mesh,'global r===')
00425       rho_prev = rho_curr
00426       ! check
00427       res_norm = Vect_dot_product(r,r)
00428       !write (stream,*) "Norm is", res_norm
00429       ratio_norm = res_norm / init_norm
00430 
00431       time_iterations = time_iterations + MPI_WTime()-t1
00432 
00433       if (ismaster()) &
00434            write(stream, '(i5,a,e22.15)') it,': res_norm=',dsqrt(res_norm)
00435     end do
00436 
00437     if(pstream/=0) then
00438        write(pstream, "(I0,':pcg iterations:',I0)") myrank, it
00439        write(pstream, "(I0,':iterations time:',F0.3)") myrank, time_iterations
00440        write(pstream, "(I0,':preconditioner time:',F0.3)") myrank, time_preconditioner
00441     end if
00442 
00443     if (ismaster()) then
00444       call CalculateEigenvalues(it,dd,ee,alpha,beta,maxit)
00445       cond_num=dd(it)/dd(1)
00446       write(stream,'(a,i3,a,e10.4)') '#it:',it,' Cond#: ',cond_num
00447     endif
00448     !deallocate(beta)
00449     !deallocate(alpha)
00450     !deallocate(ee)
00451     !deallocate(dd)
00452     ! Deallocate auxiliary data structures
00453     ! helped to assist with pmvm
00454  !  call pmvmCommStructs_destroy()
00455 
00456     !deallocate(b_k)
00457     !deallocate(z)
00458     !deallocate(m)
00459     !deallocate(q)
00460     !deallocate(p)
00461     !deallocate(r)
00462   end subroutine pcg_weigs
00463 
00464 
00465   subroutine CalculateEigenvalues(it,dd,ee,alpha,beta,MaxIt)
00466     ! Finds eigenvalues using Lanczos connection
00467     implicit none
00468     integer :: it,MaxIt
00469     real(kind=rk),dimension(:),pointer :: dd,ee,alpha,beta
00470     integer :: i,j,ierr
00471 
00472     dd(1) = 1.0_rk/alpha(1)
00473     ee(1) = 0.0_rk
00474     do i=2,it
00475       if (alpha(i)==0.or.beta(i)<0) then
00476         print *,'Unable to compute eigenvalue number',i,' and '
00477         do j=1,i
00478           print *,j,' alpha:',alpha(j),' beta:',beta
00479         enddo
00480         return
00481       else
00482         dd(i) = 1.0_rk/alpha(i) + beta(i)/alpha(i-1)
00483         ee(i) = -dsqrt(beta(i))/alpha(i-1)
00484       endif
00485     enddo
00486     !call tql1(it,dd,ee,MaxIt,ierr)
00487     call tql1(it,dd,ee,ierr)
00488     if (ierr > 0) then
00489       print *,'Cancelled after EigMaxIt while calculating eig',ierr
00490     endif
00491   end subroutine CalculateEigenvalues
00492 
00493   subroutine tql1(n,d,e,ierr)
00494     implicit none
00495     integer :: i,j,l,m,n,ii,l1,l2,mml,ierr
00496     real(kind=rk) :: d(n),e(n)
00497     real(kind=rk) :: c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2 !,pythag
00498 !
00499 !       this subroutine is a translation of the algol procedure tql1,
00500 !       num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
00501 !       wilkinson.
00502 !       handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
00503 !
00504 !       this subroutine finds the eigenvalues of a symmetric
00505 !       tridiagonal matrix by the ql method.
00506 !
00507 !       on input
00508 !          n is the order of the matrix.
00509 !          d contains the diagonal elements of the input matrix.
00510 !          e contains the subdiagonal elements of the input matrix
00511 !            in its last n-1 positions.  e(1) is arbitrary.
00512 !
00513 !        on output
00514 !          d contains the eigenvalues in ascending order.  if an
00515 !            error exit is made, the eigenvalues are correct and
00516 !            ordered for indices 1,2,...ierr-1, but may not be
00517 !            the smallest eigenvalues.
00518 !          e has been destroyed.
00519 !          ierr is set to
00520 !            zero       for normal return,
00521 !            j          if the j-th eigenvalue has not been
00522 !                       determined after 30 iterations.
00523 !       calls pythag for  sqrt(a*a + b*b) .
00524 !
00525 !       questions and comments should be directed to burton s. garbow,
00526 !       mathematics and computer science div, argonne national laboratory
00527 !
00528 !       this version dated august 1983.
00529 !
00530 !       ------------------------------------------------------------------
00531 !
00532         ierr = 0
00533         if (n .eq. 1) go to 1001
00534 !
00535         do 100 i = 2, n
00536     100 e(i-1) = e(i)
00537 !
00538         f = 0.0e0
00539         tst1 = 0.0e0
00540         e(n) = 0.0e0
00541 
00542         do 290 l = 1, n
00543            j = 0
00544            h = abs(d(l)) + abs(e(l))
00545            if (tst1 .lt. h) tst1 = h
00546 !       .......... look for small sub-diagonal element ..........
00547            do 110 m = l, n
00548               tst2 = tst1 + abs(e(m))
00549               if (tst2 .eq. tst1) go to 120
00550 !       .......... e(n) is always zero, so there is no exit
00551 !                  throu     gh the bottom of the loop ..........
00552     110    continue
00553 
00554     120    if (m .eq. l)      go to 210
00555     130    if (j .eq. 30     ) go to 1000
00556            j = j + 1
00557 !       .......... form shift ..........
00558            l1 = l + 1
00559            l2 = l1 + 1
00560            g = d(l)
00561            p = (d(l1) - g) / (2.0e0 * e(l))
00562            r = pythag(p,1.0e0_rk)
00563            d(l) = e(l) / (p + sign(r,p))
00564            d(l1) = e(l) * (p + sign(r,p))
00565            dl1 = d(l1)
00566            h = g - d(l)
00567            if (l2 .gt. n) go to 145
00568 
00569            do 140 i = l2, n
00570     140    d(i) = d(i) - h
00571 
00572     145    f = f + h
00573 !       .......... ql transformation ..........
00574            p = d(m)
00575            c = 1.0e0
00576            c2 = c
00577            el1 = e(l1)
00578            s = 0.0e0
00579            mml = m - l
00580 !       .......... for i=m-1 step -1 until l do -- ..........
00581            do 200 ii = 1, mml
00582               c3 = c2
00583               c2 = c
00584               s2 = s
00585               i = m - ii
00586               g = c * e(i)
00587               h = c * p
00588               r = pythag(p,e(i))
00589               e(i+1) = s * r
00590               s = e(i) / r
00591               c = p / r
00592               p = c * d(i) - s * g
00593               d(i+1) = h + s * (c * g + s * d(i))
00594     200    continue
00595 
00596            p = -s * s2 * c3 * el1 * e(l) / dl1
00597            e(l) = s * p
00598            d(l) = c * p
00599            tst2 = tst1 + abs(e(l))
00600            if (tst2 .gt. tst1) go to 130
00601     210    p = d(l) + f
00602 !       .......... order eigenvalues ..........
00603            if (l .eq. 1) go to 250
00604 !       .......... for i=l step -1 until 2 do -- ..........
00605            do 230 ii = 2, l
00606               i = l + 2 - ii
00607               if (p .ge. d(i-1)) go to 270
00608               d(i) = d(i-1)
00609     230    continue
00610 
00611     250    i = 1
00612     270    d(i) = p
00613     290 continue
00614 
00615         go to 1001
00616 !       .......... set error -- no convergence to an
00617 !                  eigenvalue after 30 iterations ..........
00618    1000 ierr = l
00619    1001 return
00620   end subroutine tql1
00621 
00622   function pythag(a,b) result(pyt)
00623     ! finds dsqrt(a**2+b**2) without overflow or destructive underflow
00624     implicit none
00625     real(kind=rk) :: a,b,pyt
00626     real(kind=rk) :: p,r,s,t,u
00627 
00628     p = dmax1(dabs(a),dabs(b))
00629     if (p .eq. 0.0d0) go to 20
00630     r = (dmin1(dabs(a),dabs(b))/p)**2
00631  10 continue
00632     t = 4.0d0 + r
00633     if (t .eq. 4.0d0) go to 20
00634     s = r/t
00635     u = 1.0d0 + 2.0d0*s
00636     p = u*p
00637     r = (s/u)**2 * r
00638     go to 10
00639  20 pyt = p
00640     return
00641   end function pythag
00642 
00643 
00644 end module pcg_mod
00645 
00646