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 )
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
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
430 subroutine aprod(n,x,y)
432 integer(ip),
intent(in) :: n
433 real(
dp),
intent(in) :: x(n)
434 real(
dp),
intent(out) :: y(n)
437 subroutine msolve(n,x,y)
439 integer(ip),
intent(in) :: n
440 real(
dp),
intent(in) :: x(n)
441 real(
dp),
intent(out) :: y(n)
442 end subroutine msolve
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_
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
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 =
' '
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. ', &
483 'beta2 = 0. If M = I, b and x are eigenvectors of A. ', &
484 'beta1 = 0. The exact solution is x = 0. ', &
485 'A solution to (poss. singular) Ax = b found, given rtol. ', &
486 'A solution to (poss. singular) Ax = b found, given eps. ', &
487 'Pseudoinverse solution for singular LS problem, given rtol. ', &
488 'Pseudoinverse solution for singular LS problem, given eps. ', &
489 'The iteration limit was reached. ', &
490 'The operator defined by Aprod appears to be unsymmetric. ', &
491 'The operator defined by Msolve appears to be unsymmetric. ', &
492 'The operator defined by Msolve appears to be indefinite. ', &
493 'xnorm has exceeded maxxnorm or will exceed it next iteration. ', &
494 'Acond has exceeded Acondlim or 0.1/eps. ', &
495 'Least-squares problem but no converged solution yet. ', &
496 'A null vector obtained, given rtol. ' /)
498 character(len=*),
parameter :: ddebugstr1 =
"(a, T5, i0, a, 5(e12.3))"
499 character(len=*),
parameter :: ddebugstr2 =
"(5(a, i0, a, e12.3, a))"
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 )"
526 if (
present(shift))
then
534 if (
present(disable))
then
540 if (
present(itnlim))
then
547 filename =
"MINRESQLP_tmp.txt"
550 if (
present(nout))
then
552 inquire(unit=nout, opened=connected, named=named_file, name=filename)
554 if (.not. connected)
then
555 write(*,*)
"File unit 'nout' is not open."
556 if (nout==5 .or. nout == 6)
then
562 if (.not. connected)
then
563 write(*,*)
'nout_ = ', nout_
564 open(nout_, file=filename, status=
'unknown', iostat=ios)
565 write(*,*)
'ios = ', ios
567 write(*,*)
"Error opening file '", filename,
"'."
572 if (
present(rtol))
then
581 if (
present(maxxnorm))
then
582 maxxnorm_ = min(maxxnorm, one/eps)
587 if (
present(trancond))
then
588 trancond_ = min(trancond, normmax)
593 if (
present(acondlim))
then
594 acondlim_ = min(acondlim, epsinv)
599 if (
present(msolve))
then
613 beta1 =
dnrm2(n, b, 1)
630 write(nout_, headerstr) enter, n, beta1, precon_, itnlim_, rtol_, &
631 shift_, maxxnorm_, acondlim_, trancond_
642 call msolve( n, b, y )
645 beta1 = dot_product(b, y)
647 if (beta1 < zero .and.
dnrm2(n, y, 1) > eps)
then
651 if (beta1 == zero)
then
655 beta1 = sqrt( beta1 )
658 write(*,ddebugstr1)
' y_', itn_,
' = ', (y(j), j=1,nprint)
659 write(*,*)
'beta1 ', beta1
665 if (checka_ .and. precon_)
then
666 call msolve( n, y, r2 )
667 s = dot_product( y, y )
668 t = dot_product(r1, r2)
670 epsa = (abs(s) + eps) * eps**0.33333_dp
680 call aprod ( n, y, w )
681 call aprod ( n, w, r2 )
682 s = dot_product( w, w )
683 t = dot_product( y, r2)
685 epsa = (abs(s) + eps) * eps**0.33333_dp
723 headlines = 30 * lines
725 relres = rnorm_ / (anorm_*xnorm_ + beta1)
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)
755 do while (istop_ <= flag0)
776 call aprod ( n, v, y )
777 if (shift_ /= zero)
then
781 y = y + (- beta/betal) * r1
783 alfa = dot_product(v, y)
784 y = y + (- alfa/beta) * r2
788 if ( .not. precon_ )
then
789 betan =
dnrm2(n, y, 1)
791 call msolve( n, r2, y )
792 betan = dot_product(r2, y)
793 if (betan > zero)
then
795 elseif (
dnrm2(n, y, 1) > eps )
then
804 pnorm =
dnrm2(2, vec2, 1)
809 pnorm =
dnrm2(3, vec3, 1)
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
831 dlta = cs * dbar + sn * alfa
833 gbar = sn * dbar - cs * alfa
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
853 call symortho(gbar, betan, cs, sn, gama)
859 axnorm = sqrt( axnorm**2 + tau**2 );
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
877 dlta_tmp = sr2 * vepln - cr2 * dlta
878 veplnl = cr2 * vepln + sr2 * dlta
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
898 call symortho(gamal, dlta, cr1, sr1, gamal_tmp)
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
924 ul2 = ( taul2 - etal2 * ul4 - veplnl2 * ul3 ) / gamal2
926 write(*,ddebugstr2)
'tau_', itn_-2,
' = ', taul2,
', ', &
927 'eta_', itn_-2,
' = ', etal2,
', ', &
928 'vepln_', itn_-2,
' = ', veplnl2,
', ', &
929 'gama_', itn_-2,
' = ', gamal2
933 ul = ( taul - etal * ul3 - veplnl * ul2) / gamal
935 write(*,ddebugstr2)
'tau_', itn_-1,
' = ', taul,
', ', &
936 'eta_', itn_-1,
' = ', etal,
', ', &
937 'vepln_', itn_-1,
' = ', veplnl,
', ', &
938 'gamal_', itn_-1,
' = ', gamal
945 xnorm_tmp =
dnrm2(3, vec3, 1)
946 if (abs(gama) > eps)
then
948 write(*,ddebugstr2)
'tau_', itn_,
' = ', tau,
', ', &
949 'eta_', itn_,
' = ', eta,
', ', &
950 'vepln_', itn_,
' = ', vepln,
', ', &
951 'gama_', itn_,
' = ', gama
953 u = (tau - eta*ul2 - vepln*ul) / gama
954 likels = relaresl < relresl
957 if (likels .and.
dnrm2(2, vec2, 1) > maxxnorm_)
then
967 xl2norm =
dnrm2(2, vec2, 1)
971 xnorm_ =
dnrm2(3, vec3, 1)
973 if (acond_ < trancond_ .and. istop_ == flag0 .and. qlpiter == 0)
then
976 if (gama_tmp > eps)
then
978 w = (v - epln*wl2 - dlta_qlp*wl) * s
981 if (xnorm_ < maxxnorm_)
then
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)
1000 qlpiter = qlpiter + 1;
1001 if (qlpiter == 1)
then
1005 wl2 = gamal3*wl2 + veplnl2*wl + etal*w
1008 wl = gamal_qlp*wl + vepln_qlp*w
1011 xl2 = x - ul_qlp*wl - u_qlp*w
1019 else if (itn_ == 2)
then
1027 wl2 = cr2*wl2 + sr2*v
1034 x = xl2 + ul *wl + u*w
1038 write(*,*)
'MINRESQLP updates'
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_
1061 write(*,*)
'Compute the next right reflection P{', itn_-1, itn_+1,
'}'
1062 write(*,ddebugstr2)
'gama_', itn_-1,
' = ', gamal,
', ', &
1063 'epln_', itn_+1,
' = ', eplnn
1067 call symortho(gamal_tmp, eplnn, cr2, sr2, gamal)
1070 write(*,ddebugstr2)
'cr2_', itn_+1,
' = ', cr2,
', ', &
1071 'sr2_', itn_+1,
' = ', sr2,
', ', &
1072 'gama_', itn_-1,
' = ', gamal
1077 gamal_qlp = gamal_tmp
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
1097 abs_gama = abs(gama)
1099 anorm_ = max(anorm_, pnorm, gamal, abs_gama)
1103 else if (itn_ > 1)
then
1109 gmin = min(gminl2, gamal, abs_gama)
1112 acond_ = anorm_ / gmin
1115 if (istop_ /= 14) rnorm_ = phi
1116 relres = rnorm_ / (anorm_ * xnorm_ + beta1)
1119 rootl =
dnrm2(2, vec2, 1)
1120 arnorml = rnorml * rootl
1121 relaresl = rootl / anorm_
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_
1137 epsx = anorm_*xnorm_*eps
1138 if (istop_ == flag0 .or. istop_ == 14)
then
1142 if (t1 <= one )
then
1144 else if (t2 <= one )
then
1146 else if (relres <= rtol_ )
then
1148 else if (relaresl <= rtol_ )
then
1150 else if (epsx >= beta1 )
then
1152 else if (xnorm_ >= maxxnorm_)
then
1154 else if (acond_ >= acondlim_ .or. acond_ >= ieps)
then
1156 else if (itn_ >= itnlim_ )
then
1158 else if (betan < eps )
then
1162 if (disable_ .and. itn_ < itnlim_)
then
1165 if (axnorm < rtol_*anorm_*xnorm_)
then
1171 if (istop_ /= flag0)
then
1173 if (istop_ == 6 .or. istop_ == 7 .or. istop_ == 12 .or. istop_ == 13)
then
1183 call aprod ( n, x, r1 )
1184 r1 = b - r1 + shift_*x
1185 call aprod ( n, r1, wl2 )
1186 wl2 = wl2 - shift_*r1
1187 arnorm_ =
dnrm2(n, wl2, 1)
1188 if (rnorm_ > zero .and. anorm_ > zero)
then
1189 relares = arnorm_ / (anorm_*rnorm_)
1193 if (nout_ > 0 .and. .not. lastiter .and. mod(itn_-1,lines) == 0)
then
1194 if (itn_ == 101)
then
1196 headlines = 30*lines
1197 else if (itn_ == 1001)
then
1199 headlines = 30*lines
1202 if (qlpiter == 1)
then
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.
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)")
1231 write(*,*)
'istop = ', istop_
1233 if (istop_ /= flag0)
exit
1241 if (
present(istop))
then
1245 if (
present(itn))
then
1249 if (
present(rnorm))
then
1253 if (
present(arnorm))
then
1257 if (
present(xnorm))
then
1261 if (
present(anorm))
then
1265 if (
present(acond))
then
1270 write(nout_, itnstr) itn_, x(1), xnorm_, &
1271 rnorm_, arnorm_, relres, relares, anorm_, acond_, qlpstr
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_)
1336 real(dp),
intent(in) :: a, b
1337 real(dp),
intent(out) :: c, s, r
1342 logical,
parameter :: debug = .false.
1344 real(dp) :: abs_a, abs_b
1349 if (abs_b <= realmin)
then
1358 else if (abs_a <= realmin)
then
1363 else if (abs_b > abs_a)
then
1365 s = (b / abs_b) / sqrt(one + t**2)
1371 c = (a / abs_a) / sqrt(one + t**2)
1377 write(*,*)
'c = ', c,
', s = ', s,
', r = ', r
real(dp) function, public dnrm2(n, x, incx)
DNRM2 returns the euclidean norm of a vector.
Defines precision and range in real(kind=dp) and integer(kind=ip) for portability and a few constants...
real(dp), parameter, public eps
real(dp), parameter, public one
integer, parameter, public prcsn
real(dp), parameter, public zero
integer, parameter, public ip
integer, parameter, public dp
real(dp), parameter, public realmin
MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite a...
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.
subroutine, public symortho(a, b, c, s, r)
SymOrtho: Stable Householder reflection.