trlcore.f90

用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵

! $Id: trlcore.f90,v 1.1 2004/05/04 06:29:35 fowlkes Exp $
!!!
!\Description:
! the actual work routine of restarted Lanczos program for real
! symmetric eigenvalue problems
!
! user may directly invoke this sunroutine but she/he is responsible
! for input correct arguments as needed
!
! 1) info needs to be initialized
! 2) if info%nec>0, the first nec elements of eval and first nec
!    columns of evec must contain valid eigenpairs
! 3) workspace of currect size
!    eval(mev)
!    evec(lde, mev) (lde >= nrow, mev >= ned)
!    base(ldb, info%maxlan-mev+1) (ldb>=nrow, not used if mev>maxlan)
!    wrk(lwrk) minimum amount of memory required by TRLANCZOS is
!    maxlan*(maxlan+10)
! 4) if log files are to be written, the user need to open files on IO
!    unit log_io so that the log gile may be written correctly.
! 5) the user must set the timing variable info%clk_tot and
!    info%clk_max using system_clock function call in order for this
!    subroutine to track timing results correctly
!
! the operator that defines the eigenvalue problem is expected to have
! the following interface
!  Subroutine OP(nrow, ncol, xin, ldx, yout, ldy)
!    Integer, Intent(in) :: nrow, ncol, ldx, ldy
!    Real(8), Dimension(ldx*ncol), Intent(in) :: xin
!    Real(8), Dimension(ldy*ncol), Intent(out) :: yout
!  End Subroutine OP
!
!\Algorithm
!  0. initialize input vector
!  1. do while (more eigenvalues to compute .and. more MATVEC allowed)
!  2.    first step
!     o   alpha(k+1) = dot_product(q_{k+1}, Aq_{k+1})
!     o   rr = A*q_{k+1}-alpha(k+1)*q_{k+1}-\sum_{i=1}^{k} beta(i)q_i
!     o   (re-orthogonalization)
!  3.    do j = k+2, m
!     o     rr = Aq_j
!     o     alpha(j) = dot_product(q_j, rr)
!     o     rr = rr - alpha(j)*q_j - beta(j-1)*q_{j-1}
!     o     (re-orthogonalization)
!        end do j = k+2, m
!  4.    restarting
!     o   call dstqrb to decompose the tridiagonal matrix
!     o   perform convergence test
!     o   determine what and how many Ritz pairs to save
!     o   compute the Ritz pairs to be saved
!     end do while
!
! The re-orthogonalization procedure is implemented in trl_orth.  it
! produces a normalized vector rr that is guaranteed to be orthogonal
! to the current basis.  An error will be reported if it can not
! achieve its goal.
!!!
Subroutine trlanczos(op, info, nrow, mev, eval, evec, lde, base, ldb,&
     & nbas, wrk, lwrk)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nrow, ldb, lde, mev, nbas, lwrk
  Real(8) :: eval(mev)
  Real(8), Target :: evec(lde,mev), base(1:ldb,nbas), wrk(lwrk)
  Interface
     Subroutine op(nrow, ncol, xx, ldx, yy, ldy)
       Implicit None
       Integer, Intent(in) :: nrow, ncol, ldx, ldy
       Real(8), Intent(in) :: xx(ldx*ncol)
       Real(8), Intent(out) :: yy(ldy*ncol)
     End Subroutine op
     ! this routine normalizes rr and make sure it is orthogonal to all
     ! colums of v1 and v2
     Subroutine trl_orth(nrow, v1, ld1, m1, v2, ld2, m2, rr, kept,&
          & alpha, beta, wrk, lwrk, info)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T) :: info
       Integer, Intent(in) :: nrow, ld1, ld2, m1, m2, kept, lwrk
       Real(8), Intent(in), Target :: v1(ld1,m1), v2(ld2,m2)
       Real(8) :: rr(nrow), alpha(m1+m2), beta(m1+m2), wrk(lwrk)
     End Subroutine trl_orth
     Subroutine trl_check_orth(info, nrow, v1, ld1, j1, v2, ld2, j2, wrk, lwrk)
       Use trl_info
       Implicit None
       Integer, Intent(in) :: nrow, ld1, ld2, j1, j2, lwrk
       Type(TRL_INFO_T), Intent(in) :: info
       Real(8), Intent(in) :: v1(ld1,j1), v2(ld2,j2)
       Real(8) :: wrk(lwrk)
     End Subroutine trl_check_orth
     ! compute the errors in the Lanczos recurrence
     Subroutine trl_check_recurrence(op, info, nrow, v1, ld1, m1, v2,&
          & ld2, m2, kept, alpha, beta, wrk, lwrk)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nrow, ld1, ld2, m1, m2, kept, lwrk
       Real(8), Intent(in), Target :: v1(ld1,m1), v2(ld2,m2)
       Real(8), Intent(in) :: alpha(m1+m2-1), beta(m1+m2-1)
       Real(8), Target :: wrk(lwrk)
       Interface
          Subroutine op(nrow, ncol, xx, ldx, yy, ldy)
            Implicit None
            Integer, Intent(in) :: nrow, ncol, ldx, ldy
            Real(8), Intent(in) :: xx(ldx*ncol)
            Real(8), Intent(out) :: yy(ldy*ncol)
          End Subroutine op
       End Interface
     End Subroutine trl_check_recurrence
     Subroutine trl_initial_guess(nrow, evec, lde, mev, base, ldb, nbas,&
          & alpha, beta, j1, j2, info, wrk, lwrk)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T) :: info
       Integer, Intent(in) :: nrow, lde, ldb, mev, nbas, lwrk
       Integer, Intent(out) :: j1, j2
       Real(8) :: wrk(lwrk), base(ldb,nbas), evec(lde,mev)
       Real(8) :: alpha(mev+nbas-1), beta(mev+nbas-1)
     End Subroutine trl_initial_guess
     Subroutine trl_g_sum(mpicom, n, xx, wrk)
       Implicit None
       Integer, Intent(in) :: mpicom, n
       Real(8), Intent(out) :: wrk(n)
       Real(8), Intent(inout) :: xx(n)
     End Subroutine trl_g_sum
     Subroutine trl_tridiag(nd, alpha, beta, rot, alfrot, betrot, wrk,&
          & lwrk, ierr)
       Implicit None
       Integer, Intent(in) :: nd, lwrk
       Integer, Intent(out) :: ierr
       Real(8), Intent(in), Dimension(nd) :: alpha, beta
       Real(8) :: rot(nd*nd), alfrot(nd), betrot(nd), wrk(lwrk)
     End Subroutine trl_tridiag
     Subroutine trl_get_eval(nd, locked, alpha, beta, lambda, res, wrk,&
          & lwrk, ierr)
       Implicit None
       Integer, Intent(in) :: nd, locked, lwrk
       Integer :: ierr
       Real(8), Intent(in) :: alpha(nd), beta(nd)
       Real(8), Intent(out) :: lambda(nd), res(nd), wrk(lwrk)
     End Subroutine trl_get_eval
     Subroutine trl_convergence_test(nd, lambda, res, info, wrk)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T) :: info
       Integer, Intent(in) :: nd
       Real(8), Intent(in) :: lambda(nd), res(nd)
       Real(8), Intent(out) :: wrk(nd+nd)
     End Subroutine trl_convergence_test
     Subroutine trl_print_progress(info)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
     End Subroutine trl_print_progress
     Subroutine trl_shuffle_eig(nd, mnd, lambda, res, info, kept)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, mnd
       Integer :: kept
       Real(8) :: lambda(nd), res(nd)
     End Subroutine trl_shuffle_eig
     Subroutine trl_get_tvec(nd, alpha, beta, irot, nrot, rot, nlam,&
          & lambda, yy, iwrk, wrk, lwrk, ierr)
       Implicit None
       Integer, Intent(in) :: nd, nrot, nlam, irot, lwrk
       Integer :: ierr, iwrk(nd,4)
       Real(8), Intent(in) :: alpha(nd), beta(nd), lambda(nlam), rot(nrot*nrot)
       Real(8) :: wrk(lwrk), yy(nd*nlam)
     End Subroutine trl_get_tvec
     Subroutine trl_get_tvec_a(nd, kept, alpha, beta, nlam, lambda,&
          & yy, wrk, lwrk, iwrk, ierr)
       Implicit None
       Integer, Intent(in) :: kept, nd, nlam, lwrk
       Integer :: ierr, iwrk(nd)
       Real(8), Intent(in) :: lambda(nd)
       Real(8) :: alpha(nd), beta(nd), yy(nd*nd), wrk(lwrk)
     End Subroutine trl_get_tvec_a
     Subroutine trl_ritz_vectors(nrow, lck, ny, yy, ldy, vec1, ld1, m1,&
          & vec2, ld2, m2, wrk, lwrk)
       Implicit None
       Integer, Intent(in) :: nrow, ny, lck, ldy, ld1, m1, m2, ld2, lwrk
       Real(8), Intent(in) :: yy(ldy, ny)
       Real(8) :: vec1(ld1,m1), vec2(ld2,m2), wrk(lwrk)
     End Subroutine trl_ritz_vectors
  End Interface
  !
  ! local variables
  !
  Character, Parameter :: notrans = 'N'
  Real(8), Parameter :: zero=0.0D0, one=1.0D0
  Character(132) :: title
  Integer :: i1, i2, j1, j2, jnd, jml, j1n, j2n, kept, prek
  Integer :: next_test, clk1, lwrk2, chkpnt, locked
  Integer, Dimension(:), Pointer :: iwrk
  Real(8), Dimension(:), Pointer :: alpha, beta, rr, rot,&
       & alfrot, betrot, lambda, res, yy, qa, qb, wrk2
  External dgemv
  !
  ! alpha, beta, alfrot and betrot have fixed locations in wrk
  !
  Nullify(alpha, beta, rr, rot, alfrot, betrot, lambda, res, yy, qa,&
       & qb, wrk2)
  title = ''
  clk1 = 0
  i1 = info%maxlan + 1
  i2 = info%maxlan + info%maxlan
  alpha => wrk(1:info%maxlan)
  beta  => wrk(i1:i2)
  i1 = i2 + 1
  i2 = i2 + info%maxlan
  alfrot => wrk(i1:i2)
  i1 = i2 + 1
  i2 = i2 + info%maxlan
  betrot => wrk(i1:i2)
  Allocate(iwrk(info%maxlan*4))
  iwrk = 0
  If (info%cpflag.Le.0) Then
     chkpnt = info%maxmv + info%maxlan
  Else
     chkpnt = info%maxmv / info%cpflag
  End If
  locked = info%nec
  !
  ! assign values to alpha, beta
  ! uses the assumption that the content of eval(1:nec) are eigenvalues
  ! and their residual norms are zero
  !
  locked = info%nec
  If (locked.Gt.0) Then
     alpha(1:locked) = eval(1:locked)
     beta(1:locked) = zero
  End If
  ! get valid initial guess for the Lanczos iterations
  wrk2 => wrk(i2+1:lwrk)
  lwrk2 = lwrk - i2
  Call trl_initial_guess(nrow, evec, lde, mev, base, ldb, nbas, alpha,&
       & beta, j1, j2, info, wrk2, lwrk2)
  If (info%stat .Ne. 0) Goto 888
  jnd = j1 + j2
  kept = jnd
  ! we will perform the first convergence test after next_test
  ! matrix-vector multiplications
  i1 = info%ned - jnd
  next_test = i1 + Min(i1, 6, info%ned/2)
  !!*************************************************!!
  !! the restart loop --                             !!
  !! restart if there is more eigenvalues to compute !!
  !! and more MATVEC allowed to do it                !!
  !!*************************************************!!
  Do While (info%matvec.Lt.info%maxmv .And. info%nec.Lt.info%ned)
     jnd = jnd + 1
     i2 = lwrk - (jnd-locked)**2
     rot => wrk(i2+1:lwrk)
     i1 = 4*info%maxlan + 1
     wrk2 => wrk(i1:i2)              ! workspace for orthogonalization
     lwrk2 = i2 - i1 + 1
     i1 = Max(5*info%maxlan-3*locked, 4*info%maxlan)
     If (lwrk2 .Lt. i1) Then
        info%stat = -11
        Return
     End If
     !
     ! first step
     !
     If (j1 .Lt. mev) Then
        j1 = j1 + 1
        qa => evec(1:nrow, j1)
     Else
        j2 = j2+1
        qa => base(1:nrow, j2)
     End If
     If (j1 .Lt. mev) Then
        j1n = j1 + 1
        j2n = 0
     Else
        j1n = mev
        j2n = j2 + 1
     End If
     Nullify(qb)
     If (j1.Lt.mev) Then
        rr => evec(1:nrow, j1n)
     Else
        rr => base(1:nrow, j2n)
     End If
     ! perform matrix-vector multiplication
     ! record the total time and the time in MATVEC
     Call System_clock(clk1)
     If (clk1 .Le. info%clk_tot) Then
        info%tick_t = info%tick_t + (info%clk_max-info%clk_tot)
        info%tick_t = info%tick_t + clk1
        info%clk_tot = clk1
     End If
!!$print *, 'input to OP #', info%matvec + 1
!!$print *, qa
!!$print *, 'should be the same as'
!!$if (j2 .eq. 0) then
!!$print *, evec(1:nrow,j1)
!!$else
!!$print *, base(1:nrow,j2)
!!$end if
     Call OP(nrow, 1, qa, nrow, rr, nrow)
!!$print *, 'return of OP #', info%matvec + 1
!!$print *, rr
     Call add_clock_ticks(info%clk_op, info%tick_o)
!!$     Call hpcheck(info%stat)
!!$     Print *, 'hpcheck (trlcore.f90:298) returns ', info%stat,&
!!$          & ', nloop=', info%nloop, ', nec=', info%nec, ', jnd=',jnd
     info%matvec = info%matvec+1
     ! perform the Lanczos orthogonalization
     alpha(jnd) = Dot_product(qa, rr)
     Call trl_g_sum(info%mpicom, 1, alpha(jnd:jnd), wrk2)
     beta(jnd) = alpha(jnd)
     info%flop = info%flop + nrow + nrow
     If (j1 .Gt. 2) Then
        Call dgemv(notrans, nrow, j1, -one, evec, lde,&
             & beta, 1, one, rr, 1)
        info%flop = info%flop + 2*j1*nrow
     Else If (j1 .Eq. 1) Then
        rr = rr - alpha(1)*qa
        info%flop = info%flop + nrow + nrow
     Else If (j1 .Eq. 2) Then
        rr = rr - beta(1)*evec(1:nrow,1) - beta(2)*evec(1:nrow,2)
        info%flop = info%flop + 4*nrow
     End If
     If (j2 .Gt. 2) Then
        Call dgemv(notrans, nrow, j2, -one, base, ldb,&
             & beta(j1+1), 1, one, rr, 1)
        info%flop = info%flop + 2*j2*nrow
     Else If (j2 .Eq. 1) Then
        rr = rr - beta(jnd)*qa
        info%flop = info%flop + nrow + nrow
     Else If (j2 .Eq. 2) Then
        rr = rr - beta(j1+1)*base(1:nrow,1) - beta(jnd)*base(1:nrow,2)
        info%flop = info%flop + 4*nrow
     End If
     ! perform re-orthogonalization
     info%flop_h = info%flop_h - info%flop
     Call System_clock(clk1)
     Call trl_orth(nrow, evec, lde, j1, base, ldb, j2, rr,&
          & kept, alpha, beta, wrk2, lwrk2, info)
     If (info%verbose.Gt.8 .Or. info%stat.Ne.0) Then
        ! check orthogonality after the initilization step
        Call trl_check_orth(info, nrow, evec, lde, j1n, base, ldb, j2n,&
             & wrk2, lwrk2)
     End If
     Call add_clock_ticks(info%clk_orth, info%tick_h)
     info%flop_h = info%flop_h + info%flop
     If (info%stat .Ne. 0) Goto 888
     If (info%verbose .Gt. 5) Call print_alpha_beta
     !
     ! transform the matrix formed by alpha and beta into a
     ! tridiagonal matrix, rot stores the transformation matrix
     !
     alfrot(1:locked) = alpha(1:locked)
     betrot(1:locked) = zero
     i1 = jnd-locked
     Call trl_tridiag(i1, alpha(locked+1:locked+i1), beta(locked+1:locked+i1),&
          & rot, alfrot(locked+1:locked+i1), betrot(locked+1:locked+i1), wrk2,&
          & lwrk2, info%stat)
     info%flop = info%flop + 8*i1*i1*i1/3  ! Golub:1996:MC, P415
     If (info%stat .Ne. 0) Goto 888
     betrot(jnd) = beta(jnd)
     !!******************************************************!!
     !! regular iterations of Lanczos algorithm (inner loop) !!
     !!******************************************************!!
     Do While (jnd .Lt. info%klan .And. info%nec .Lt. info%ned)
        qb => qa
        qa => rr
        j1 = j1n
        j2 = j2n
        jnd = jnd + 1
        If (j1n .Lt. mev) Then
           j1n = j1n + 1
        Else
           j2n = j2n + 1
        End If
        If (j1.Lt.mev) Then
           rr => evec(1:nrow, j1n)
        Else
           rr => base(1:nrow, j2n)
        End If
        ! MATVEC
        Call System_clock(clk1)
        Call OP(nrow, 1, qa, nrow, rr, nrow)
        Call add_clock_ticks(info%clk_op, info%tick_o)