Millepede-II  V04-08-03
minresModule.f90
Go to the documentation of this file.
1 
3 
4 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5 ! File minresModule.f90
6 !
38 
40 
41  use minresdatamodule, only : dp
42 
43  implicit none
44  public :: minres
45 
46 contains
47 
48  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
49 
294  subroutine minres( n, Aprod, Msolve, b, shift, checkA, precon, &
295  x, itnlim, nout, rtol, &
296  istop, itn, Anorm, Acond, rnorm, Arnorm, ynorm )
297 
298  integer, intent(in) :: n, itnlim, nout
299  logical, intent(in) :: checka, precon
300  real(dp), intent(in) :: b(n)
301  real(dp), intent(in) :: shift, rtol
302  real(dp), intent(out) :: x(n)
303  integer, intent(out) :: istop, itn
304  real(dp), intent(out) :: anorm, acond, rnorm, arnorm, ynorm
305 
306  interface
307  subroutine aprod (n,x,y) ! y := A*x
308  use minresdatamodule
309  integer, intent(in) :: n
310  real(dp), intent(in) :: x(n)
311  real(dp), intent(out) :: y(n)
312  end subroutine aprod
313 
314  subroutine msolve(n,x,y) ! Solve M*y = x
315  use minresdatamodule
316  integer, intent(in) :: n
317  real(dp), intent(in) :: x(n)
318  real(dp), intent(out) :: y(n)
319  end subroutine msolve
320  end interface
321 
322 ! Local arrays and variables
323  real(dp) :: r1(n), r2(n), v(n), w(n), w1(n), w2(n), y(n)
324  real(dp) :: alfa , beta , beta1 , cs , &
325  dbar , delta , denom , diag , &
326  eps , epsa , epsln , epsr , epsx , &
327  gamma , gbar , gmax , gmin , &
328  oldb , oldeps, qrnorm, phi , phibar, &
329  rhs1 , rhs2 , rnorml, rootl , &
330  arnorml, relarnorml, &
331  s , sn , t , tnorm2, ynorm2, z
332  integer :: i
333  logical :: debug, prnt
334 
335  ! Local constants
336  real(dp), parameter :: zero = 0.0, one = 1.0
337  real(dp), parameter :: ten = 10.0
338  character(len=*), parameter :: enter = ' Enter MINRES. '
339  character(len=*), parameter :: exitt = ' Exit MINRES. '
340  character(len=*), parameter :: msg(-1:8) = &
341  (/ 'beta2 = 0. If M = I, b and x are eigenvectors of A', & ! -1
342  'beta1 = 0. The exact solution is x = 0 ', & ! 0
343  'Requested accuracy achieved, as determined by rtol ', & ! 1
344  'Reasonable accuracy achieved, given eps ', & ! 2
345  'x has converged to an eigenvector ', & ! 3
346  'Acond has exceeded 0.1/eps ', & ! 4
347  'The iteration limit was reached ', & ! 5
348  'Aprod does not define a symmetric matrix ', & ! 6
349  'Msolve does not define a symmetric matrix ', & ! 7
350  'Msolve does not define a pos-def preconditioner ' /) ! 8
351  !-------------------------------------------------------------------
352 
353  intrinsic :: abs, dot_product, epsilon, min, max, sqrt
354 
355  ! Print heading and initialize.
356 
357  debug = .false.
358  eps = epsilon(eps)
359  if (nout > 0) then
360  write(nout, 1000) enter, n, checka, precon, itnlim, rtol, shift
361  end if
362  istop = 0
363  itn = 0
364  anorm = zero
365  acond = zero
366  rnorm = zero
367  ynorm = zero
368  x(1:n) = zero
369  arnorml = zero
370  gmin = zero
371  gmax = zero
372 
373  !-------------------------------------------------------------------
374  ! Set up y and v for the first Lanczos vector v1.
375  ! y = beta1 P' v1, where P = C**(-1).
376  ! v is really P' v1.
377  !-------------------------------------------------------------------
378  y = b
379  r1 = b
380  if ( precon ) call msolve( n, b, y )
381  beta1 = dot_product(b,y)
382 
383  if (beta1 < zero) then ! M must be indefinite.
384  istop = 8
385  go to 900
386  end if
387 
388  if (beta1 == zero) then ! b = 0 exactly. Stop with x = 0.
389  istop = 0
390  go to 900
391  end if
392 
393  beta1 = sqrt( beta1 ) ! Normalize y to get v1 later.
394 
395  !-------------------------------------------------------------------
396  ! See if Msolve is symmetric.
397  !-------------------------------------------------------------------
398  if (checka .and. precon) then
399  call msolve( n, y, r2 )
400  s = dot_product(y ,y )
401  t = dot_product(r1,r2)
402  z = abs(s - t)
403  epsa = (s + eps) * eps**0.33333
404  if (z > epsa) then
405  istop = 7
406  go to 900
407  end if
408  end if
409 
410  !-------------------------------------------------------------------
411  ! See if Aprod is symmetric. Initialize Arnorm.
412  !-------------------------------------------------------------------
413  if (checka) then
414  call aprod ( n, y, w )
415  call aprod ( n, w, r2 )
416  s = dot_product(w,w )
417  t = dot_product(y,r2)
418  z = abs(s - t)
419  epsa = (s + eps) * eps**0.33333
420  if (z > epsa) then
421  istop = 6
422  go to 900
423  end if
424  arnorml = sqrt(s);
425  else
426  call aprod ( n, y, w )
427  arnorml = sqrt( dot_product(w,w) )
428  end if
429 
430  !-------------------------------------------------------------------
431  ! Initialize other quantities.
432  !-------------------------------------------------------------------
433  oldb = zero
434  beta = beta1
435  dbar = zero
436  epsln = zero
437  qrnorm = beta1
438  phibar = beta1
439  rhs1 = beta1
440  rhs2 = zero
441  tnorm2 = zero
442  ynorm2 = zero
443  cs = - one
444  sn = zero
445  w(1:n) = zero
446  w2(1:n)= zero
447  r2(1:n)= r1
448 
449  if (debug) then
450  write(*,*) ' '
451  write(*,*) 'b ', b
452  write(*,*) 'beta ', beta
453  write(*,*) ' '
454  end if
455 
456  !===================================================================
457  ! Main iteration loop.
458  !===================================================================
459  do
460  itn = itn + 1 ! k = itn = 1 first time through
461 
462  !----------------------------------------------------------------
463  ! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
464  ! The general iteration is similar to the case k = 1 with v0 = 0:
465  !
466  ! p1 = Operator * v1 - beta1 * v0,
467  ! alpha1 = v1'p1,
468  ! q2 = p2 - alpha1 * v1,
469  ! beta2^2 = q2'q2,
470  ! v2 = (1/beta2) q2.
471  !
472  ! Again, y = betak P vk, where P = C**(-1).
473  ! .... more description needed.
474  !----------------------------------------------------------------
475  s = one / beta ! Normalize previous vector (in y).
476  v = s*y(1:n) ! v = vk if P = I
477 
478  call aprod ( n, v, y )
479  y = y - shift*v ! call daxpy ( n, (- shift), v, 1, y, 1 )
480  if (itn >= 2) then
481  y = y - (beta/oldb)*r1 ! call daxpy ( n, (- beta/oldb), r1, 1, y, 1 )
482  end if
483 
484  alfa = dot_product(v,y) ! alphak
485  y = y - (alfa/beta)*r2 ! call daxpy ( n, (- alfa/beta), r2, 1, y, 1 )
486  r1 = r2
487  r2 = y
488  if ( precon ) call msolve( n, r2, y )
489 
490  oldb = beta ! oldb = betak
491  beta = dot_product(r2,y) ! beta = betak+1^2
492  if (beta < zero) then
493  istop = 6
494  go to 900
495  end if
496 
497  beta = sqrt( beta ) ! beta = betak+1
498  tnorm2 = tnorm2 + alfa**2 + oldb**2 + beta**2
499 
500  if (itn == 1) then ! Initialize a few things.
501  if (beta/beta1 <= ten*eps) then ! beta2 = 0 or ~ 0.
502  istop = -1 ! Terminate later.
503  end if
504  !tnorm2 = alfa**2
505  gmax = abs( alfa ) ! alpha1
506  gmin = gmax ! alpha1
507  end if
508 
509  ! Apply previous rotation Qk-1 to get
510  ! [deltak epslnk+1] = [cs sn][dbark 0 ]
511  ! [gbar k dbar k+1] [sn -cs][alfak betak+1].
512 
513  oldeps = epsln
514  delta = cs * dbar + sn * alfa ! delta1 = 0 deltak
515  gbar = sn * dbar - cs * alfa ! gbar 1 = alfa1 gbar k
516  epsln = sn * beta ! epsln2 = 0 epslnk+1
517  dbar = - cs * beta ! dbar 2 = beta2 dbar k+1
518 
519  ! Compute the next plane rotation Qk
520 
521  gamma = sqrt( gbar**2 + beta**2 ) ! gammak
522  cs = gbar / gamma ! ck
523  sn = beta / gamma ! sk
524  phi = cs * phibar ! phik
525  phibar = sn * phibar ! phibark+1
526 
527  if (debug) then
528  write(*,*) ' '
529  write(*,*) 'v ', v
530  write(*,*) 'alfa ', alfa
531  write(*,*) 'beta ', beta
532  write(*,*) 'gamma', gamma
533  write(*,*) 'delta', delta
534  write(*,*) 'gbar ', gbar
535  write(*,*) 'epsln', epsln
536  write(*,*) 'dbar ', dbar
537  write(*,*) 'phi ', phi
538  write(*,*) 'phiba', phibar
539  write(*,*) ' '
540  end if
541 
542  ! Update x.
543 
544  denom = one/gamma
545 
546  do i = 1, n
547  w1(i) = w2(i)
548  w2(i) = w(i)
549  w(i) = ( v(i) - oldeps*w1(i) - delta*w2(i) ) * denom
550  x(i) = x(i) + phi * w(i)
551  end do
552 
553  ! Go round again.
554 
555  gmax = max( gmax, gamma )
556  gmin = min( gmin, gamma )
557  z = rhs1 / gamma
558  ynorm2 = z**2 + ynorm2
559  rhs1 = rhs2 - delta * z
560  rhs2 = - epsln * z
561 
562  ! Estimate various norms and test for convergence.
563 
564  anorm = sqrt( tnorm2 )
565  ynorm = sqrt( ynorm2 )
566  epsa = anorm * eps
567  epsx = anorm * ynorm * eps
568  epsr = anorm * ynorm * rtol
569  diag = gbar
570  if (diag == zero) diag = epsa
571 
572  qrnorm = phibar
573  rnorml = rnorm
574  rnorm = qrnorm
575  rootl = sqrt( gbar**2 +dbar**2 ) ! norm([gbar; dbar]);
576  arnorml = rnorml*rootl ! ||A r_{k-1} ||
577  relarnorml = rootl / anorm; ! ||Ar|| / (||A|| ||r||)
578  !relArnorml = Arnorml / Anorm; ! ||Ar|| / ||A||
579 
580  ! Estimate cond(A).
581  ! In this version we look at the diagonals of R in the
582  ! factorization of the lower Hessenberg matrix, Q * H = R,
583  ! where H is the tridiagonal matrix from Lanczos with one
584  ! extra row, beta(k+1) e_k^T.
585 
586  acond = gmax / gmin
587 
588  ! See if any of the stopping criteria are satisfied.
589  ! In rare cases, istop is already -1 from above (Abar = const*I).
590 
591  if (istop == 0) then
592  if (itn >= itnlim ) istop = 5
593  if (acond >= 0.1d+0/eps) istop = 4
594  if (epsx >= beta1 ) istop = 3
595  ! original
596  !if (qrnorm <= epsx .or. relArnorml <= epsx) istop = 2
597  !if (qrnorm <= epsr .or. relArnorml <= epsr) istop = 1
598  ! C. Kleinwort, DESY, 131002
599  if (qrnorm <= epsx .or. relarnorml <= eps) istop = 2
600  if (qrnorm <= epsr .or. relarnorml <= rtol) istop = 1
601  end if
602 
603 
604  ! See if it is time to print something.
605 
606  if (nout > 0) then
607  prnt = .false.
608  if (n <= 40 ) prnt = .true.
609  if (itn <= 10 ) prnt = .true.
610  if (itn >= itnlim - 10) prnt = .true.
611  if (mod(itn,10) == 0) prnt = .true.
612  if (qrnorm <= ten * epsx) prnt = .true.
613  if (qrnorm <= ten * epsr) prnt = .true.
614  if (relarnorml<= ten*epsx) prnt = .true.
615  if (relarnorml<= ten*epsr) prnt = .true.
616  if (acond >= 1.0d-2/eps ) prnt = .true.
617  if (istop /= 0 ) prnt = .true.
618 
619  if ( prnt ) then
620  if ( itn == 1) write(nout, 1200)
621  write(nout, 1300) itn, x(1), qrnorm, anorm, acond
622  if (mod(itn,10) == 0) write(nout, 1500)
623  end if
624  end if
625  if (istop /= 0) exit
626 
627  end do
628  !===================================================================
629  ! End of iteration loop.
630  !===================================================================
631 
632  ! Display final status.
633 
634 900 arnorm = arnorml
635  if (nout > 0) then
636  write(nout, 2000) exitt, istop, itn, &
637  exitt, anorm, acond, &
638  exitt, rnorm, ynorm, arnorm
639  write(nout, 3000) exitt, msg(istop)
640  end if
641 
642  return
643 
644  1000 format(// 1p, a, 5x, 'Solution of symmetric Ax = b' &
645  / ' n =', i7, 5x, 'checkA =', l4, 12x, &
646  'precon =', l4 &
647  / ' itnlim =', i7, 5x, 'rtol =', e11.2, 5x, &
648  'shift =', e23.14)
649  1200 format(// 5x, 'itn', 8x, 'x(1)', 10x, &
650  'norm(r)', 3x, 'norm(A)', 3x, 'cond(A)')
651  1300 format(1p, i8, e19.10, 3e10.2)
652  1500 format(1x)
653  2000 format(/ 1p, a, 5x, 'istop =', i3, 14x, 'itn =', i8 &
654  / a, 5x, 'Anorm =', e12.4, 5x, 'Acond =', e12.4 &
655  / a, 5x, 'rnorm =', e12.4, 5x, 'ynorm =', e12.4, 5x, 'Arnorml =', e12.4)
656  3000 format( a, 5x, a )
657 
658  end subroutine minres
659 
660 end module minresmodule
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
Definition: mpnum.f90:2443
Defines real(kind=dp) and a few constants for use in other modules.
integer, parameter, public dp
MINRES solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite and/...
subroutine, public minres(n, Aprod, Msolve, b, shift, checkA, precon, x, itnlim, nout, rtol, istop, itn, Anorm, Acond, rnorm, Arnorm, ynorm)
Solution of linear equation system.