Millepede-II  V04-07-04
minresqlpModule.f90
Go to the documentation of this file.
1 
3 
4 
5 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 ! File minresqlpModule.f90
7 !
40 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41 
43 
45  use minresqlpblasmodule, only : dnrm2
46 
47  implicit none
48 
49  private ! sets default for module
50  public :: minresqlp, symortho
51 
52 contains
53 
54  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
55 
410 
411  subroutine minresqlp( n, Aprod, b, shift, Msolve, disable, nout, &
412  itnlim, rtol, maxxnorm, trancond, Acondlim, &
413  x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond )
414  ! Inputs
415  integer(ip), intent(in) :: n
416  real(dp), intent(in) :: b(n)
417  integer(ip), intent(in), optional :: itnlim, nout
418  logical, intent(in), optional :: disable
419  real(dp), intent(in), optional :: shift
420  real(dp), intent(in), optional :: rtol, maxxnorm, trancond, acondlim
421 
422  ! Outputs
423  real(dp), intent(out) :: x(n)
424  integer(ip), intent(out), optional :: istop, itn
425  real(dp), intent(out), optional :: rnorm, arnorm, xnorm, anorm, acond
426 
427  optional :: msolve
428 
429  interface
430  subroutine aprod(n,x,y) ! y := A*x
432  integer(ip), intent(in) :: n
433  real(dp), intent(in) :: x(n)
434  real(dp), intent(out) :: y(n)
435  end subroutine aprod
436 
437  subroutine msolve(n,x,y) ! Solve M*y = x
439  integer(ip), intent(in) :: n
440  real(dp), intent(in) :: x(n)
441  real(dp), intent(out) :: y(n)
442  end subroutine msolve
443  end interface
444 
445  intrinsic :: abs, sqrt, present, floor, log10, dot_product
446 
447  ! Local arrays and variables
448  real(dp) :: shift_
449  real(dp) :: rtol_, maxxnorm_, trancond_, acondlim_
450  real(dp) :: rnorm_, arnorm_, xnorm_, anorm_, acond_
451  logical :: checka_, precon_, disable_
452  integer(ip) :: itnlim_, nout_, istop_, itn_
453 
454  real(dp) :: r1(n), r2(n), v(n), w(n), wl(n), wl2(n),&
455  xl2(n), y(n), vec2(2), vec3(3)
456  real(dp) :: abs_gama, acondl, alfa, anorml, arnorml,&
457  axnorm, beta, beta1, betal, betan, &
458  epsa, epsx, gminl2, ieps, pnorm, &
459  relares, relaresl, relres, relresl, &
460  rnorml, rootl, t1, t2, xl2norm, &
461  xnorm_tmp, xnorml, z, &
462  cr1, cr2, cs, dbar, dlta, dlta_qlp, &
463  dlta_tmp, dltan, epln, eplnn, eta, etal,&
464  etal2, gama, gama_qlp, gama_tmp, gamal, &
465  gamal_qlp, gamal_tmp, gamal2, gamal3, &
466  gbar, gmin, gminl, phi, s, sn, sr1, sr2,&
467  t, tau, taul, taul2, u, u_qlp, &
468  ul, ul_qlp, ul2, ul2_qlp, ul3, ul4, &
469  vepln, vepln_qlp, veplnl, veplnl2, x1last
470 
471  integer(ip) :: j, qlpiter, headlines, lines, nprint, flag0, ios
472  logical :: prnt, done, lastiter, connected, named_file, likels
473  character(len=20) :: filename
474  character(len=2) :: qlpstr = ' '
475 
476  ! Local constants
477  real(dp), parameter :: epsinv = 10.0_dp**floor(log10(one/eps))
478  real(dp) :: normmax = 10.0_dp**floor(log10(one/eps)/2)
479  character(len=*), parameter :: enter = ' Enter MINRES-QLP. '
480  character(len=*), parameter :: exitt = ' Exit MINRES-QLP. '
481  character(len=*), parameter :: msg(1:15) = &
482  (/ 'beta_{k+1} < eps. ', & ! 1
483  'beta2 = 0. If M = I, b and x are eigenvectors of A. ', & ! 2
484  'beta1 = 0. The exact solution is x = 0. ', & ! 3
485  'A solution to (poss. singular) Ax = b found, given rtol. ', & ! 4
486  'A solution to (poss. singular) Ax = b found, given eps. ', & ! 5
487  'Pseudoinverse solution for singular LS problem, given rtol. ', & ! 6
488  'Pseudoinverse solution for singular LS problem, given eps. ', & ! 7
489  'The iteration limit was reached. ', & ! 8
490  'The operator defined by Aprod appears to be unsymmetric. ', & ! 9
491  'The operator defined by Msolve appears to be unsymmetric. ', & ! 10
492  'The operator defined by Msolve appears to be indefinite. ', & ! 11
493  'xnorm has exceeded maxxnorm or will exceed it next iteration. ', & ! 12
494  'Acond has exceeded Acondlim or 0.1/eps. ', & ! 13
495  'Least-squares problem but no converged solution yet. ', & ! 14
496  'A null vector obtained, given rtol. ' /) ! 15
497 
498  character(len=*), parameter :: ddebugstr1 = "(a, T5, i0, a, 5(e12.3))"
499  character(len=*), parameter :: ddebugstr2 = "(5(a, i0, a, e12.3, a))"
500 
501  character(len=*), parameter :: headerstr = &
502  "(// 1p, a, 4x, 'Solution of symmetric Ax = b'" // &
503  " / ' n =', i7, 6x, '||b|| =', e11.2, 3x," // &
504  " 'precon =', l4 " // &
505  " / ' itnlim =', i7, 6x, 'rtol =', e11.2, 3x," // &
506  " 'shift =', e23.14 " // &
507  " / ' maxxnorm =', e11.2, 2x, 'Acondlim =', e11.2, 3x,"// &
508  " 'trancond =', e11.2)"
509  character(len=*), parameter :: tableheaderstr = &
510  "(// ' iter x(1) xnorm rnorm Arnorm '," // &
511  " 'Compatible LS norm(A) cond(A)')"
512  character(len=*), parameter :: itnstr = "(1p, i8, e19.10, 7e10.2, a)"
513  character(len=*), parameter :: finalstr1 = &
514  "(/ 1p, a, 5x, a, i3, 14x, a, i8 " // &
515  " / a, 5x, a, e12.4, 5x, a, e12.4 " // &
516  " / a, 5x, a, e12.4, 5x, a, e12.4 " // &
517  " / a, 5x, a, e12.4 )"
518  character(len=*), parameter :: finalstr2 = "( a, 5x, a )"
519 
520  !------------------------------------------------------------------
521  prnt = .false.
522  t1 = zero
523  t2 = zero
524  gminl = zero
525  ! Optional inputs
526  if (present(shift)) then
527  shift_ = shift
528  else
529  shift_ = zero
530  end if
531 
532  checka_ = .true.
533 
534  if (present(disable)) then
535  disable_ = disable
536  else
537  disable_ = .false.
538  end if
539 
540  if (present(itnlim)) then
541  itnlim_ = itnlim
542  else
543  itnlim_ = 4 * n
544  end if
545 
546  connected = .false.
547  filename = "MINRESQLP_tmp.txt"
548  nout_ = 10
549 
550  if (present(nout)) then
551  nout_ = nout
552  inquire(unit=nout, opened=connected, named=named_file, name=filename)
553  !write(*,*) connected, named_file, filename
554  if (.not. connected) then
555  write(*,*) "File unit 'nout' is not open."
556  if (nout==5 .or. nout == 6) then
557  nout_ = 10
558  end if
559  end if
560  end if
561 
562  if (.not. connected) then
563  write(*,*) 'nout_ = ', nout_
564  open(nout_, file=filename, status='unknown', iostat=ios)
565  write(*,*) 'ios = ', ios
566  if (ios /= 0) then
567  write(*,*) "Error opening file '", filename, "'."
568  stop
569  end if
570  end if
571 
572  if (present(rtol)) then
573  rtol_ = rtol
574  else
575  rtol_ = eps
576  end if
577 
578  if (prcsn == 6) then
579  normmax = 1.0e4_dp
580  end if
581  if (present(maxxnorm)) then
582  maxxnorm_ = min(maxxnorm, one/eps)
583  else
584  maxxnorm_ = normmax
585  end if
586 
587  if (present(trancond)) then
588  trancond_ = min(trancond, normmax)
589  else
590  trancond_ = normmax
591  end if
592 
593  if (present(acondlim)) then
594  acondlim_ = min(acondlim, epsinv)
595  else
596  acondlim_ = epsinv
597  end if
598 
599  if (present(msolve)) then
600  precon_ = .true.
601  else
602  precon_ = .false.
603  end if
604 
605  !------------------------------------------------------------------
606  ! Print heading and initialize.
607  !------------------------------------------------------------------
608  nprint = min(n,20)
609  !debug = .true.
610  lastiter = .false.
611  flag0 = 0
612  istop_ = flag0
613  beta1 = dnrm2(n, b, 1)
614  ieps = 0.1_dp/eps
615  itn_ = 0
616  qlpiter = 0
617  xnorm_ = zero
618  xl2norm = zero
619  axnorm = zero
620  anorm_ = zero
621  acond_ = one
622  pnorm = zero
623  relresl = zero
624  relaresl = zero
625  x = zero
626  xl2 = zero
627  x1last = x(1)
628 
629  if (nout_ > 0) then
630  write(nout_, headerstr) enter, n, beta1, precon_, itnlim_, rtol_, &
631  shift_, maxxnorm_, acondlim_, trancond_
632  end if
633 
634  !------------------------------------------------------------------
635  ! Set up y and v for the first Lanczos vector v1.
636  ! y = beta1 P'v1, where P = C**(-1).
637  ! v is really P'v1.
638  !------------------------------------------------------------------
639  y = b
640  r1 = b
641  if ( precon_ ) then
642  call msolve( n, b, y )
643  end if
644 
645  beta1 = dot_product(b, y)
646 
647  if (beta1 < zero .and. dnrm2(n, y, 1) > eps) then ! M must be indefinite.
648  istop_ = 11
649  end if
650 
651  if (beta1 == zero) then ! b = 0 exactly. Stop with x = 0.
652  istop_ = 3
653  end if
654 
655  beta1 = sqrt( beta1 ) ! Normalize y to get v1 later.
656 
657  if (debug) then
658  write(*,ddebugstr1) ' y_', itn_, ' = ', (y(j), j=1,nprint)
659  write(*,*) 'beta1 ', beta1
660  end if
661 
662  !------------------------------------------------------------------
663  ! See if Msolve is symmetric.
664  !------------------------------------------------------------------
665  if (checka_ .and. precon_) then
666  call msolve( n, y, r2 )
667  s = dot_product( y, y )
668  t = dot_product(r1, r2)
669  z = abs( s - t )
670  epsa = (abs(s) + eps) * eps**0.33333_dp
671  if (z > epsa) then
672  istop_ = 10
673  end if
674  end if
675 
676  !------------------------------------------------------------------
677  ! See if Aprod is symmetric.
678  !------------------------------------------------------------------
679  if (checka_) then
680  call aprod ( n, y, w ) ! w = A*y
681  call aprod ( n, w, r2 ) ! r2 = A*w
682  s = dot_product( w, w )
683  t = dot_product( y, r2)
684  z = abs( s - t )
685  epsa = (abs(s) + eps) * eps**0.33333_dp
686  if (z > epsa) then
687  istop_ = 9
688  end if
689  end if
690 
691  !------------------------------------------------------------------
692  ! Initialize other quantities.
693  !------------------------------------------------------------------
694  tau = zero
695  taul = zero
696  gmin = zero
697  beta = zero
698  betan = beta1
699  dbar = zero
700  dltan = zero
701  eta = zero
702  etal = zero
703  etal2 = zero
704  vepln = zero
705  veplnl = zero
706  veplnl2= zero
707  eplnn = zero
708  gama = zero
709  gamal = zero
710  gamal2 = zero
711  phi = beta1
712  cr1 = -one
713  sr1 = zero
714  cr2 = -one
715  sr2 = zero
716  cs = -one
717  sn = zero
718  ul3 = zero
719  ul2 = zero
720  ul = zero
721  u = zero
722  lines = 1
723  headlines = 30 * lines
724  rnorm_ = betan
725  relres = rnorm_ / (anorm_*xnorm_ + beta1)
726  relares= zero
727  r2 = b
728  w = zero ! vector of zeros
729  wl = zero
730  done = .false.
731 
732  if (debug) then
733  write(*,*)
734  write(*,*) 'Checking variable values before main loop'
735  write(*,*) 'istop ', istop_, ' done ', done, ' itn ', itn_, ' QLPiter' , qlpiter
736  write(*,*) 'lastiter', lastiter, ' lines ', lines, ' headlines ', headlines
737  write(*,*) 'beta ', beta, ' tau ', tau, ' taul ', taul, ' phi ', phi
738  write(*,*) 'betan ', betan, ' gmin ', gmin, ' cs ', cs, ' sn ', sn
739  write(*,*) 'cr1 ', cr1, ' sr1 ', sr1, ' cr2 ', cr2, ' sr2 ', sr2
740  write(*,*) 'dltan ', dltan, ' eplnn ', eplnn, ' gama ', gama, ' gamal ', gamal
741  write(*,*) 'gamal2 ', gamal2, ' eta ', eta, ' etal ', etal, 'etal2', etal2
742  write(*,*) 'vepln ', vepln, ' veplnl', veplnl, ' veplnl2 ', veplnl2, ' ul3 ', ul3
743  write(*,*) 'ul2 ', ul2, ' ul ', ul, ' u ', u, ' rnorm ', rnorm_
744  write(*,*) 'xnorm ', xnorm_, ' xl2norm ', xl2norm, ' pnorm ', pnorm, ' Axnorm ', axnorm
745  write(*,*) 'Anorm ', anorm_, ' Acond ', acond_, ' relres ', relres
746  write(*,ddebugstr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
747  write(*,ddebugstr1) 'w_', itn_, ' = ', (w(j), j=1,nprint)
748  write(*,ddebugstr1) 'x_', itn_, ' = ', (x(j), j=1,nprint)
749  end if
750 
751 
752  !------------------------------------------------------------------
753  ! Main iteration loop.
754  !------------------------------------------------------------------
755  do while (istop_ <= flag0)
756  itn_ = itn_ + 1 ! k = itn = 1 first time through
757 
758  !-----------------------------------------------------------------
759  ! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
760  ! The general iteration is similar to the case k = 1 with v0 = 0:
761  !
762  ! p1 = Operator * v1 - beta1 * v0,
763  ! alpha1 = v1'p1,
764  ! q2 = p2 - alpha1 * v1,
765  ! beta2^2 = q2'q2,
766  ! v2 = (1/beta2) q2.
767  !
768  ! Again, y = betak P vk, where P = C**(-1).
769  ! .... more description needed.
770  !-----------------------------------------------------------------
771 
772  betal = beta; ! betal = betak
773  beta = betan;
774  s = one / beta ! Normalize previous vector (in y).
775  v = s*y; ! v = vk if P = I.
776  call aprod ( n, v, y )
777  if (shift_ /= zero) then
778  y = y - shift_ * v
779  end if
780  if (itn_ >= 2) then
781  y = y + (- beta/betal) * r1
782  end if
783  alfa = dot_product(v, y) ! alphak
784  y = y + (- alfa/beta) * r2
785  r1 = r2
786  r2 = y
787 
788  if ( .not. precon_ ) then
789  betan = dnrm2(n, y, 1) ! betan = ||y||_2
790  else
791  call msolve( n, r2, y )
792  betan = dot_product(r2, y) ! betan = betak+1^2
793  if (betan > zero) then
794  betan = sqrt(betan)
795  elseif ( dnrm2(n, y, 1) > eps ) then ! M must be indefinite.
796  istop_ = 11
797  exit
798  end if
799  end if
800 
801  if (itn_ == 1) then
802  vec2(1) = alfa
803  vec2(2) = betan
804  pnorm = dnrm2(2, vec2, 1)
805  else
806  vec3(1) = beta
807  vec3(2) = alfa
808  vec3(3) = betan
809  pnorm = dnrm2(3, vec3, 1)
810  end if
811 
812 
813  if (debug) then
814  write(*,*)
815  write(*,*) 'Lanczos iteration ', itn_
816  write(*,ddebugstr1) 'v_', itn_, ' = ', (v(j), j=1,nprint)
817  write(*,ddebugstr1) 'r1_', itn_, ' = ', (r1(j), j=1,nprint)
818  write(*,ddebugstr1) 'r2_', itn_, ' = ', (r2(j), j=1,nprint)
819  write(*,ddebugstr1) 'y_', itn_, ' = ', (y(j), j=1,nprint)
820  write(*,ddebugstr2) 'alpha_', itn_, ' = ', alfa, ', ', &
821  'beta_', itn_, ' = ', beta, ', ', &
822  'beta_', itn_+1, ' = ', betan, ', ', &
823  'pnorm_', itn_, ' = ', pnorm
824  end if
825 
826  ! Apply previous left reflection Qk-1 to get
827  ! [deltak epslnk+1] = [cs sn][dbark 0 ]
828  ! [gbar k dbar k+1] [sn -cs][alfak betak+1].
829 
830  dbar = dltan
831  dlta = cs * dbar + sn * alfa ! dlta1 = 0 deltak
832  epln = eplnn;
833  gbar = sn * dbar - cs * alfa ! gbar 1 = alfa1 gbar k
834  eplnn = sn * betan ! eplnn2 = 0 epslnk+1
835  dltan = - cs * betan ! dbar 2 = beta2 dbar k+1
836  dlta_qlp = dlta;
837 
838  if (debug) then
839  write(*,*)
840  write(*,*) 'Apply previous left reflection Q_{', itn_-1, ',', itn_, '}'
841  write(*,ddebugstr2) 'c_', itn_-1, ' = ', cs, ', ', &
842  's_', itn_-1, ' = ', sn
843  write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
844  'gbar_', itn_, ' = ', gbar, ', ', &
845  'epln_', itn_+1, ' = ', eplnn, ', ', &
846  'dbar_', itn_+1, ' = ', dltan
847  end if
848 
849  ! Compute the current left reflection Qk
850  gamal3 = gamal2
851  gamal2 = gamal
852  gamal = gama
853  call symortho(gbar, betan, cs, sn, gama)
854  gama_tmp = gama;
855  taul2 = taul
856  taul = tau
857  tau = cs * phi
858  phi = sn * phi ! phik
859  axnorm = sqrt( axnorm**2 + tau**2 );
860 
861  if (debug) then
862  write(*,*)
863  write(*,*) 'Compute the current left reflection Q_{', itn_, ',', itn_+1, '}'
864  write(*,ddebugstr2) 'c_', itn_, ' = ', cs, ', ', &
865  's_', itn_, ' = ', sn
866  write(*,ddebugstr2) 'tau_', itn_, ' = ', tau, ', ', &
867  'phi_', itn_, ' = ', phi, ', ', &
868  'gama_', itn_, ' = ', gama
869  end if
870 
871  ! Apply the previous right reflection P{k-2,k}
872 
873  if (itn_ > 2) then
874  veplnl2 = veplnl
875  etal2 = etal
876  etal = eta
877  dlta_tmp = sr2 * vepln - cr2 * dlta
878  veplnl = cr2 * vepln + sr2 * dlta
879  dlta = dlta_tmp
880  eta = sr2 * gama
881  gama = -cr2 * gama
882  if (debug) then
883  write(*,*)
884  write(*,*) 'Apply the previous right reflection P_{', itn_-2, ',', itn_, '}'
885  write(*,ddebugstr2) 'cr2_', itn_, ' = ', cr2, ', ', &
886  'sr2_', itn_, ' = ', sr2
887  write(*,ddebugstr2) 'gamal2_', itn_, ' = ', gamal2, ', ', &
888  'gamal_', itn_, ' = ', gamal, ', ', &
889  'gama_', itn_, ' = ', gama
890  write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
891  'vepln_', itn_-1, ' = ', veplnl
892  end if
893  end if
894 
895 
896  ! Compute the current right reflection P{k-1,k}, P_12, P_23,...
897  if (itn_ > 1) then
898  call symortho(gamal, dlta, cr1, sr1, gamal_tmp)
899  gamal = gamal_tmp
900  vepln = sr1 * gama
901  gama = -cr1 * gama
902 
903  if (debug) then
904  write(*,*)
905  write(*,*) 'Apply the current right reflection P_{', itn_-1, ',', itn_, '}'
906  write(*,ddebugstr2) 'cr1_ ', itn_, ' = ', cr1, ', ', &
907  'sr1_' , itn_, ' = ', sr1
908  write(*,ddebugstr2) 'gama_', itn_-2, ' = ', gamal2, ', ', &
909  'gama_', itn_-1, ' = ', gamal, ', ', &
910  'gama_', itn_, ' = ', gama
911  write(*,ddebugstr2) 'dlta_', itn_, ' = ', dlta, ', ', &
912  'vepln_', itn_-1, ' = ', veplnl, ', ', &
913  'eta_', itn_, ' = ', eta
914  end if
915  end if
916 
917  ! Update xnorm
918 
919  xnorml = xnorm_
920  ul4 = ul3
921  ul3 = ul2
922 
923  if (itn_ > 2) then
924  ul2 = ( taul2 - etal2 * ul4 - veplnl2 * ul3 ) / gamal2
925  if (debug) then
926  write(*,ddebugstr2) 'tau_', itn_-2, ' = ', taul2, ', ', &
927  'eta_', itn_-2, ' = ', etal2, ', ', &
928  'vepln_', itn_-2, ' = ', veplnl2, ', ', &
929  'gama_', itn_-2, ' = ', gamal2
930  end if
931  end if
932  if (itn_ > 1) then
933  ul = ( taul - etal * ul3 - veplnl * ul2) / gamal
934  if (debug) then
935  write(*,ddebugstr2) 'tau_', itn_-1, ' = ', taul, ', ', &
936  'eta_', itn_-1, ' = ', etal, ', ', &
937  'vepln_', itn_-1, ' = ', veplnl, ', ', &
938  'gamal_', itn_-1, ' = ', gamal
939  end if
940  end if
941 
942  vec3(1) = xl2norm
943  vec3(2) = ul2
944  vec3(3) = ul
945  xnorm_tmp = dnrm2(3, vec3, 1) ! norm([xl2norm ul2 ul]);
946  if (abs(gama) > eps) then ! .and. xnorm_tmp < maxxnorm_) then
947  if (debug) then
948  write(*,ddebugstr2) 'tau_', itn_, ' = ', tau, ', ', &
949  'eta_', itn_, ' = ', eta, ', ', &
950  'vepln_', itn_, ' = ', vepln, ', ', &
951  'gama_', itn_, ' = ', gama
952  end if
953  u = (tau - eta*ul2 - vepln*ul) / gama
954  likels = relaresl < relresl
955  vec2(1) = xnorm_tmp
956  vec2(2) = u
957  if (likels .and. dnrm2(2, vec2, 1) > maxxnorm_) then
958  u = zero
959  istop_ = 12
960  end if
961  else
962  u = zero
963  istop_ = 14
964  end if
965  vec2(1) = xl2norm
966  vec2(2) = ul2
967  xl2norm = dnrm2(2, vec2, 1)
968  vec3(1) = xl2norm
969  vec3(2) = ul
970  vec3(3) = u
971  xnorm_ = dnrm2(3, vec3, 1)
972 
973  if (acond_ < trancond_ .and. istop_ == flag0 .and. qlpiter == 0) then !! MINRES updates
974  wl2 = wl
975  wl = w
976  if (gama_tmp > eps) then
977  s = one / gama_tmp
978  w = (v - epln*wl2 - dlta_qlp*wl) * s
979  end if
980 
981  if (xnorm_ < maxxnorm_) then
982  x1last = x(1)
983  x = x + tau*w
984  else
985  istop_ = 12
986  lastiter = .true.
987  end if
988 
989  if (debug) then
990  write(*,*)
991  write(*,*) 'MINRES updates'
992  write(*,ddebugstr2) 'gama_tmp_', itn_, ' = ', gama_tmp, ', ', &
993  'tau_', itn_, ' = ', tau, ', ', &
994  'epln_', itn_, ' = ', epln, ', ', &
995  'dlta_QLP_', itn_, ' = ', dlta_qlp
996  write(*,ddebugstr1) 'v_', itn_ , ' = ', (v(j), j=1,nprint)
997  write(*,ddebugstr1) 'w_', itn_ , ' = ', (w(j), j=1,nprint)
998  end if
999  else !! MINRES-QLP updates
1000  qlpiter = qlpiter + 1;
1001  if (qlpiter == 1) then
1002  xl2 = zero ! vector
1003  if (itn_ > 1) then ! construct w_{k-3}, w_{k-2}, w_{k-1}
1004  if (itn_ > 3) then
1005  wl2 = gamal3*wl2 + veplnl2*wl + etal*w
1006  end if ! w_{k-3}
1007  if (itn_ > 2) then
1008  wl = gamal_qlp*wl + vepln_qlp*w
1009  end if ! w_{k-2}
1010  w = gama_qlp*w
1011  xl2 = x - ul_qlp*wl - u_qlp*w
1012  end if
1013  end if
1014 
1015  if (itn_ == 1) then
1016  wl2 = wl
1017  wl = sr1*v
1018  w = -cr1*v
1019  else if (itn_ == 2) then
1020  wl2 = wl
1021  wl = cr1*w + sr1*v
1022  w = sr1*w - cr1*v
1023  else
1024  wl2 = wl
1025  wl = w
1026  w = sr2*wl2 - cr2*v
1027  wl2 = cr2*wl2 + sr2*v
1028  v = cr1*wl + sr1*w
1029  w = sr1*wl - cr1*w
1030  wl = v
1031  end if
1032  x1last = x(1)
1033  xl2 = xl2 + ul2*wl2
1034  x = xl2 + ul *wl + u*w
1035 
1036  if (debug) then
1037  write(*,*)
1038  write(*,*) 'MINRESQLP updates'
1039  end if
1040  end if
1041 
1042  if (debug) then
1043  write(*,*)
1044  write(*,*) 'Update u, w, x and xnorm'
1045  write(*,ddebugstr2) 'u_', itn_-2, ' = ', ul2, ', ', &
1046  'u_', itn_-1, ' = ', ul, ', ', &
1047  'u_', itn_, ' = ', u
1048  write(*,ddebugstr1) 'w_', itn_-2, ' = ', (wl2(j), j=1,nprint)
1049  write(*,ddebugstr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
1050  write(*,ddebugstr1) 'w_', itn_, ' = ', (w(j), j=1,nprint)
1051  write(*,ddebugstr1) 'x_', itn_, ' = ', (x(j), j=1,nprint)
1052  write(*,ddebugstr2) 'xnorm_', itn_-2, ' = ', xl2norm, ', ', &
1053  'xnorm_', itn_-1, ' = ', xnorml, ', ', &
1054  'xnorm_', itn_, ' = ', xnorm_
1055  end if
1056 
1057  ! Compute the next right reflection P{k-1,k+1}
1058 
1059  if (debug) then
1060  write(*,*)
1061  write(*,*) 'Compute the next right reflection P{', itn_-1, itn_+1,'}'
1062  write(*,ddebugstr2) 'gama_', itn_-1, ' = ', gamal, ', ', &
1063  'epln_', itn_+1, ' = ', eplnn
1064  end if
1065 
1066  gamal_tmp = gamal
1067  call symortho(gamal_tmp, eplnn, cr2, sr2, gamal)
1068 
1069  if (debug) then
1070  write(*,ddebugstr2) 'cr2_', itn_+1, ' = ', cr2, ', ', &
1071  'sr2_', itn_+1, ' = ', sr2, ', ', &
1072  'gama_', itn_-1, ' = ', gamal
1073  end if
1074 
1075  ! Store quantities for transfering from MINRES to MINRES-QLP
1076 
1077  gamal_qlp = gamal_tmp
1078  vepln_qlp = vepln
1079  gama_qlp = gama
1080  ul2_qlp = ul2
1081  ul_qlp = ul
1082  u_qlp = u
1083 
1084  if (debug) then
1085  write(*,*)
1086  write(*,*) 'Store quantities for transfering from MINRES to MINRES-QLP '
1087  write(*,ddebugstr2) 'gama_QLP_', itn_-1, ' = ', gamal_qlp, ', ', &
1088  'vepln_QLP_',itn_, ' = ', vepln_qlp, ', ', &
1089  'gama_QLP_', itn_, ' = ', gama_qlp
1090  write(*,ddebugstr2) 'u_QLP_', itn_-2, ' = ', ul2_qlp, ', ', &
1091  'u_QLP_', itn_-1, ' = ', ul_qlp, ', ', &
1092  'u_QLP_', itn_, ' = ', u_qlp
1093  end if
1094 
1095  ! Estimate various norms
1096 
1097  abs_gama = abs(gama)
1098  anorml = anorm_
1099  anorm_ = max(anorm_, pnorm, gamal, abs_gama)
1100  if (itn_ == 1) then
1101  gmin = gama
1102  gminl = gmin
1103  else if (itn_ > 1) then
1104  gminl2 = gminl
1105  gminl = gmin
1106  vec3(1) = gminl2
1107  vec3(2) = gamal
1108  vec3(3) = abs_gama
1109  gmin = min(gminl2, gamal, abs_gama)
1110  end if
1111  acondl = acond_
1112  acond_ = anorm_ / gmin
1113  rnorml = rnorm_
1114  relresl = relres
1115  if (istop_ /= 14) rnorm_ = phi
1116  relres = rnorm_ / (anorm_ * xnorm_ + beta1)
1117  vec2(1) = gbar
1118  vec2(2) = dltan
1119  rootl = dnrm2(2, vec2, 1)
1120  arnorml = rnorml * rootl
1121  relaresl = rootl / anorm_
1122 
1123  if (debug) then
1124  write(*,*)
1125  write(*,*) 'Estimate various norms '
1126  write(*,ddebugstr2) 'gmin_', itn_, ' = ', gmin, ', ', &
1127  'pnorm_', itn_, ' = ', pnorm, ', ', &
1128  'rnorm_', itn_, ' = ', rnorm_, ', ', &
1129  'Arnorm_', itn_-1, ' = ', arnorml
1130  write(*,ddebugstr2) 'Axnorm_', itn_, ' = ', axnorm, ', ', &
1131  'Anorm_', itn_, ' = ', anorm_, ', ', &
1132  'Acond_', itn_, ' = ', acond_
1133  end if
1134 
1135  ! See if any of the stopping criteria are satisfied.
1136 
1137  epsx = anorm_*xnorm_*eps
1138  if (istop_ == flag0 .or. istop_ == 14) then
1139  t1 = one + relres
1140  t2 = one + relaresl
1141  end if
1142  if (t1 <= one ) then
1143  istop_ = 5 ! Accurate Ax=b solution
1144  else if (t2 <= one ) then
1145  istop_ = 7 ! Accurate LS solution
1146  else if (relres <= rtol_ ) then
1147  istop_ = 4 ! Good enough Ax=b solution
1148  else if (relaresl <= rtol_ ) then
1149  istop_ = 6 ! Good enough LS solution
1150  else if (epsx >= beta1 ) then
1151  istop_ = 2 ! x is an eigenvector
1152  else if (xnorm_ >= maxxnorm_) then
1153  istop_ = 12 ! xnorm exceeded its limit
1154  else if (acond_ >= acondlim_ .or. acond_ >= ieps) then
1155  istop_ = 13 ! Huge Acond
1156  else if (itn_ >= itnlim_ ) then
1157  istop_ = 8 ! Too many itns
1158  else if (betan < eps ) then
1159  istop_ = 1 ! Last iteration of Lanczos, rarely happens
1160  end if
1161 
1162  if (disable_ .and. itn_ < itnlim_) then
1163  istop_ = flag0
1164  done = .false.
1165  if (axnorm < rtol_*anorm_*xnorm_) then
1166  istop_ = 15
1167  lastiter = .false.
1168  end if
1169  end if
1170 
1171  if (istop_ /= flag0) then
1172  done = .true.
1173  if (istop_ == 6 .or. istop_ == 7 .or. istop_ == 12 .or. istop_ == 13) then
1174  lastiter = .true.
1175  end if
1176  if (lastiter) then
1177  itn_ = itn_ - 1
1178  acond_ = acondl
1179  rnorm_ = rnorml
1180  relres = relresl
1181  end if
1182 
1183  call aprod ( n, x, r1 )
1184  r1 = b - r1 + shift_*x ! r1 to temporarily store residual vector
1185  call aprod ( n, r1, wl2 ) ! wl2 to temporarily store A*r1
1186  wl2 = wl2 - shift_*r1
1187  arnorm_ = dnrm2(n, wl2, 1)
1188  if (rnorm_ > zero .and. anorm_ > zero) then
1189  relares = arnorm_ / (anorm_*rnorm_)
1190  end if
1191  end if
1192 
1193  if (nout_ > 0 .and. .not. lastiter .and. mod(itn_-1,lines) == 0) then
1194  if (itn_ == 101) then
1195  lines = 10
1196  headlines = 30*lines
1197  else if (itn_ == 1001) then
1198  lines = 100
1199  headlines = 30*lines
1200  end if
1201 
1202  if (qlpiter == 1) then
1203  qlpstr = ' P'
1204  else
1205  qlpstr = ' '
1206  end if
1207  end if
1208 
1209 
1210  ! See if it is time to print something.
1211 
1212  if (nout_ > 0) then
1213  prnt = .false.
1214  if (n <= 40 ) prnt = .true.
1215  if (itn_ <= 10 ) prnt = .true.
1216  if (itn_ >= itnlim_ - 10) prnt = .true.
1217  if (mod(itn_-1,10) == 0 ) prnt = .true.
1218  if (qlpiter == 1 ) prnt = .true.
1219  if (acond_ >= 0.01_dp/eps ) prnt = .true.
1220  if (istop_ /= flag0 ) prnt = .true.
1221 
1222  if ( prnt ) then
1223  if (itn_ == 1) write(nout_, tableheaderstr)
1224  write(nout_, itnstr) itn_-1, x1last, xnorml, &
1225  rnorml, arnorml, relresl, relaresl, anorml, acondl, qlpstr
1226  if (mod(itn_,10) == 0) write(nout_, "(1x)")
1227  end if
1228  end if
1229 
1230  if (debug) then
1231  write(*,*) 'istop = ', istop_
1232  end if
1233  if (istop_ /= flag0) exit
1234  enddo
1235  !===================================================================
1236  ! End of iteration loop.
1237  !===================================================================
1238 
1239  ! Optional outputs
1240 
1241  if (present(istop)) then
1242  istop = istop_
1243  end if
1244 
1245  if (present(itn)) then
1246  itn = itn_
1247  end if
1248 
1249  if (present(rnorm)) then
1250  rnorm = rnorm_
1251  end if
1252 
1253  if (present(arnorm)) then
1254  arnorm = arnorm_
1255  end if
1256 
1257  if (present(xnorm)) then
1258  xnorm = xnorm_
1259  end if
1260 
1261  if (present(anorm)) then
1262  anorm = anorm_
1263  end if
1264 
1265  if (present(acond)) then
1266  acond = acond_
1267  end if
1268 
1269  if ( prnt ) then
1270  write(nout_, itnstr) itn_, x(1), xnorm_, &
1271  rnorm_, arnorm_, relres, relares, anorm_, acond_, qlpstr
1272  end if
1273 
1274  ! Display final status.
1275 
1276  if (nout_ > 0) then
1277  write(nout_, finalstr1) exitt, 'istop =', istop_, 'itn =', itn_, &
1278  exitt, 'Anorm =', anorm_, 'Acond =', acond_, &
1279  exitt, 'rnorm =', rnorm_, 'Arnorm =', arnorm_, &
1280  exitt, 'xnorm =', xnorm_
1281  write(nout_, finalstr2) exitt, msg(istop_)
1282  end if
1283 
1284  return
1285 end subroutine minresqlp
1286 
1287  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1288 
1289  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1332  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1333 
1334  subroutine symortho(a, b, c, s, r)
1335 
1336  real(dp), intent(in) :: a, b
1337  real(dp), intent(out) :: c, s, r
1338 
1339  intrinsic :: abs, sqrt
1340 
1341  ! Local variables
1342  logical, parameter :: debug = .false.
1343  real(dp) :: t
1344  real(dp) :: abs_a, abs_b
1345 
1346  abs_a = abs(a)
1347  abs_b = abs(b)
1348 
1349  if (abs_b <= realmin) then
1350  s = zero
1351  r = abs_a
1352  if (a == zero) then
1353  c = one
1354  else
1355  c = a / abs_a
1356  end if
1357 
1358  else if (abs_a <= realmin) then
1359  c = zero;
1360  r = abs_b
1361  s = b / abs_b
1362 
1363  else if (abs_b > abs_a) then
1364  t = a / b
1365  s = (b / abs_b) / sqrt(one + t**2)
1366  c = s * t
1367  r = b / s ! computationally better than r = a / c since |c| <= |s|
1368 
1369  else
1370  t = b / a
1371  c = (a / abs_a) / sqrt(one + t**2)
1372  s = c * t;
1373  r = a / c ! computationally better than r = b / s since |s| <= |c|
1374  end if
1375 
1376  if (debug) then
1377  write(*,*) 'c = ', c, ', s = ', s, ', r = ', r
1378  end if
1379 
1380  end subroutine symortho
1381 
1382 end module minresqlpmodule
minresqlpmodule
MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite a...
Definition: minresqlpModule.f90:42
minresqlpmodule::symortho
subroutine, public symortho(a, b, c, s, r)
SymOrtho: Stable Householder reflection.
Definition: minresqlpModule.f90:1335
minresqlpdatamodule::debug
logical, public debug
Definition: minresqlpDataModule.f90:55
minresqlpdatamodule::zero
real(dp), parameter, public zero
Definition: minresqlpDataModule.f90:49
minresqlpdatamodule::realmin
real(dp), parameter, public realmin
Definition: minresqlpDataModule.f90:50
minresqlpdatamodule::ip
integer, parameter, public ip
Definition: minresqlpDataModule.f90:47
minresqlpblasmodule
Definition: minresqlpBlasModule.f90:34
minresqlpmodule::minresqlp
subroutine, public minresqlp(n, Aprod, b, shift, Msolve, disable, nout, itnlim, rtol, maxxnorm, trancond, Acondlim, x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond)
Solution of linear equation system or least squares problem.
Definition: minresqlpModule.f90:414
minresqlpdatamodule::one
real(dp), parameter, public one
Definition: minresqlpDataModule.f90:49
minresqlpblasmodule::dnrm2
real(dp) function, public dnrm2(n, x, incx)
DNRM2 returns the euclidean norm of a vector.
Definition: minresqlpBlasModule.f90:185
minresqlpdatamodule::dp
integer, parameter, public dp
Definition: minresqlpDataModule.f90:43
minresqlpdatamodule
Defines precision and range in real(kind=dp) and integer(kind=ip) for portability and a few constants...
Definition: minresqlpDataModule.f90:35
minresqlpdatamodule::prcsn
integer, parameter, public prcsn
Definition: minresqlpDataModule.f90:52
minresqlpdatamodule::eps
real(dp), parameter, public eps
Definition: minresqlpDataModule.f90:49